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Abstract — In this paper, we have computationally examined oscillating flow (zero mean) between two 
parallel plates with a sudden change in cross section. The flow was assumed to be laminar incompressible 
with the inflow velocity uniform over the channel cross section but varying sinusoidally with time. The 
cases studied cover wide ranges of Re m ,« (from 187.5 to 2000), Va (from l to 10.66), the expansion ratio 
(1:2 and 1:4) and A t (2 and 4). Also, three different geometric cases were discussed: (a) asymmetric 
expansion/contraction; (b) symmetric expansion/contraction; and (c) symmetric blunt body. For these 
oscillating flow conditions, the fluid undergoes sudden expansion in one-half of the cycle and sudden 
contraction in the other half. The instantaneous friction factor, for some ranges of Re^ and Va, deviated 
substantially from the steady-state friction factor for the same flow parameters. A region has been 
identified (see Fig. 3) below which the flow is laminar quasi-steady. A videotape showing computer 
simulations of the oscillating flow demonstrates the usefulness of the current analyses in providing 
information on the transient hydraulic phenomena. 


NOMENCLATURE 


A r — Relative amplitude of the fluid displacement 
[see equation (II)] 

D h — Hydraulic diameter of the smaller channel 
/—Instantaneous friction factor ( = 2r w /pt/* ) 
h — Height of the smaller channel 
H — Height of the larger channel 
/ — Total channel length 
P — Pressure 

Re— Instantaneous Reynolds number 
S — Step size (see Fig. I) 

St — Strouhal number [see equation (12)] 
t — Time 

T — Time period for one cycle 
U — A* -component of velocity 
V— T-component of velocity 
Va — Valensi number [see equation (10)] 


X — Distance along the channel axis 
Y— Distance normal to the channel axis 

Greek symbols 

p — Density of the fluid 
p — Dynamic viscosity of the fluid 
o — Frequency of oscillation 

* Subscripts 

Dh — Based on the hydraulic diameter 
i — Based on the inlet condition 
m — Mean value 

max — Maximum during the cycle 
min — Minimum during the cycle 
0 — Initial condition/reference state 
ss — Steady state 
w — At the wall 


INTRODUCTION 

Several engineering applications encounter unsteady flow as well as sudden changes in the channel 
cross section. In free-piston Stirling engine applications, the flow oscillates around a zero mean 
while sudden changes in the cross section take place at the components interface. As an example, 
in the NASA SPRE (Space Power Research Engine) the flow goes through a sudden change in cross 
section at the interface between the heater and the expansion space as well as between the cooler 
and the compression space. Today, a steady-state correlation for the fluid flow and heat transfer 
are used in the design analyses of such engines. 

Typically, under these oscillating flow conditions the fluid undergoes a sudden expansion in 
one-half of the cycle and a sudden contraction in the other half. The flow reversal is caused not 
only by the flow oscillation but also by the sudden expansion. 

Since the flow in Stirling engines oscillates in a cyclic manner, the velocity and temperature 
profiles differ significantly from those obtained for steady flows [1-3]. Accordingly, the friction 
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factor and heat transfer coefficient are considerably different from those of the steady-state 
correlations. 

Several computational, experimental and analytical efforts have been conducted to examine the 
oscillating flows within a straight channel (circular pipe and parallel plates); Ibrahim et al. [1] and 
Kurzweg [2] have performed several numerical investigations. Simon and Seume (3, 4] conducted 
an experimental analysis for oscillating flow in a circular pipe. Their results show a significant 
increase in the friction factor as compared to steady flow conditions. 

A literature survey shows that several investigations have been conducted for flows with a sudden 
change in the area of the cross-section. Examples involve steady unidirectional flow over a 
backward-facing step [5-8] and flow through a sudden contraction in a channel [9-1 1]. The results 
obtained showed an enhancement in the friction factor as compared to a straight geometry under 
similar flow conditions. 

In this paper, results from a computational analysis of the flow between two parallel plates with 
a sudden change in cross section are presented. The flow parameters are selected to emulate the 
NASA SPRE. Several geometries and expansion ratios have been examined to identify the flow 
characteristics under different engine operating conditions. 

Since experimental data for such a problem do not exist, a careful step-by-step procedure 
has been used for validating the computer code. This includes, examination of false diffu- 
sion, comparison with available data for steady flows as well as solving for impulsively started 
flows. 


ANALYSIS 

Assumptions 

Figure 1 shows the channel with two parallel plates and a sudden change in cross section 
and the Cartesian coordinate system used for the present analyses. Different geometries are 
considered: (a) asymmetric expansion/contraction [Fig. 1(a)]; (b) symmetric expansion/ 
contraction [Fig. 1(b)]; and (c) a symmetric blunt body [Fig. 1(c)]. The following assumptions 
were made: (1) the flow is laminar incompressible with constant thermophysical properties; 
(2) the inlet velocity is uniform but varies sinusoidally with time; and (3) the location of 
the step is far away from both ends; this is chosen to isolate the effect of the step on the flow 
field from the end effects. 



(a) Asymmetric expansion/contraction 



(b) Symmetric expansion/contraction 



(c) Symmetric blunt body 

Fig. 1. The different geometries examined: (a) asymmetric expansion/contraction; (b) symmetric 
expansion/contraction; (c) symmetric blunt body. 
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Governing Equations 

For unsteady laminar flow the following system of equations is employed. 


Continuity equation 


d(pu) djpv ) _ 
dx dy 


Momentum equations 

Two equations for momentum result from conservation of momentum in the x- and y -directions, 
respectively: 
x -momentum, 

djpu) _ d(pu 2 + x„) dp 

dt dx dy dx * w 

and 

y -momentum, 

djpv) d(puv + x yx ) d(pv 2 + t yy ) dp 

dt dx dy dy ' W 

The momentum equations are cast into a standardized form by utilizing the relationship between 
the viscous stress and the strain rates: 


(du 8v \ 

~ Xxy ~ ~ Tyx = ti \dy + d^)‘ 

Equations of the following form result, which are parabolic in time but elliptic in space coordinates: 

aw , ‘’("“’-'‘I) /(»— »£) » | 

1 . I ^ ^ ** l 1 


where 


d(pu) 


dx dy dy 


* j i : 

dx dy 


•isl , is) 

t I 


For incompressible/constant property flows, the source terms are zero from the continuity 
equation. 

Dimensionless Parameters 

Different dimensionless parameters characterize the unsteady flow in the channel under 
consideration: 

(1) Re„„; for oscillating flows the mean flow velocity for a cycle is zero, therefore the Reynolds 
number is based on the maximum amplitude of the velocity during each cycle: 

Re_ — 


( 9 ) 
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(2) Va; the frequency of oscillation has been expressed in dimensionless form as the Valensi 
number: 


Va = 


CD(d ) 2 



( 10 ) 


(3) A x \ the relative amplitude of fluid displacement is defined as the maximum fluid displacement 
during half a cycle divided by the channel length, based on the assumption that the fluid moves 
as a slug flow through the passage: 


A r 



(ID 


Three different physical situations can be identified: 

(a) A r < 1; part of the fluid oscillates within the passage without exiting. 

(b) A r = 1; the volume of fluid displaced in half a cycle is exactly equal to the volume of fluid 
contained within the passage. 

(c) A r > 1; the volume of fluid displaced during half a cycle is greater than the volume of fluid 
contained within the passage. 


(4) St; the Strouhal number is a combination of Re^ and Va: 


St = 4.0* 



( 12 ) 


It should be noted that for the same channel geometry, a constant A v would also imply a 
constant St. 


(5) Another important physical parameter is the channel expansion ratio h ///(sometimes referred 
to in terms of the step height). 


Boundary Conditions 

The following boundary conditions are applied: 


(1) Solid walls, 

(2) Axis of symmetry, 

(3) Inlet plane 


U = V = 0 . 



u =0. 


Ui = U m sin (cot). 


(13) 

(14) 

(15) 


(4) Outlet plane; the exit plane is chosen to be sufficiently far away from the zones of recirculation, 
therefore the gradients normal to the exit plane (i.e. along the streamwise direction) can be 
neglected: 


5m ^ _ 

dx dx 


(16) 


It should be noted that for oscillating flows the inlet and outlet planes are switched at the 
appropriate time step so that a flow reversal is implemented numerically. 


Numerical Method 

The analysis performed utilizes a modified version of the computer code CAST, developed by 
Peric and Scheuerer [12]. The original code is capable of solving two-dimensional, steady and 
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Table l . The different cases studied in the present work 


Test 

case 

Expansion 

ratio 

R e™, 

Va 

No. of 
axes of 
symmetry 


r w 

Mesh 

A 

1 

2 

187.5 

1.0 

N/A 

2.0 

1.2 

84x22 

8 

\ 

2 

187.5 

1.0 

1 

4.0 

1.2 

84x22 

C 

1 

2 

187.5 

1.0 

2 

4.0 

1.2 

84 x 22 

D 

1 

4 

187.5 

1.0 

N/A 

2.0 

1.2 

84x22 

E 

1 

2 

1000.0 

5.33 

N/A 

2.0 

1.2 

84x22 

F 

1 

2 

1000.0 

5.33 

l 

4.0 

1.2 

84x22 

G 

1 

2 

1000.0 

5.33 

2 

4.0 

1.2 

84x22 

H 

1 

2 

1000.0 

5.33 

1 

2.0 

1.2 

94x22 

I 

l 

2 

1000.0 

5.33 

2 

2.0 

1.2 

94x22 

J 

1 

4 

1000.0 

5.33 

N/A 

2.0 

1.2 

84x32 

K 

I 

2 

2000.0 

10.66 

N/A 

2.0 

1.2 

84x22 

L 

1 

4 

2000.0 

10.66 

N/A 

2.0 

1.2 

84x32 


unsteady, unidirectional flow problems. It has been modified to handle time-dependent boundary 
conditions for oscillating flows. CAST solves two-dimensional Navier-Stokes equations for 
laminar flows utilizing a collocated grid. A special velocity-pressure coupling [12] is used, based 
on the staggered grid concept to prevent oscillatory pressure solution [13]. The numerical solution 
procedure is a conservative finite volume method using primitive variables such as velocities, 
pressure and enthalpy. The basic principal involved in this method is to balance the dependent 
variable fluxes at the inlet and outlet of each control volume within the analysis domain. The 
solution procedure employed is the well-known SIMPLE algorithm by Patankar [13]. 

For all cases investigated, the flow cycle was divided into 60 time steps of 6-degree intervals. At 
least 3 cycles were run for each case to achieve a converged solution. A 0.2% convergence criteria 
was used in this study. The CPU time required on a Cray X-MP/Y-MP ranged from 3600 to 8000 s 
(for 3 cycles), depending upon the size of the mesh used (see Table 1). 


CODE VALIDATION 

Several computational experiments were conducted to validate the CAST code as listed below: 

1. The code predictions for the reattachment length and the minimum and maximum velocities at 
various locations along the channel axis, for Re ~ 50 and 150 and expansion ratios of 1:1-2 and 
l : l .5, were compared with similar steady flow results by Morgan et al. [8]; the comparisons were 
good. 

2. Comparisons were made between the present code predictions and the numerical computations 
by Chiu [6] for steady flow over a backward-facing step, asymmetric channel with a 1:1.5 
expansion ratio and Re = 916. The present code prediction for the friction factor is within 5% 
of Chiu’s results. 

3. Comparison was made, for the size of the separation bubble before the step in a forward- facing 
step flow at Re = 200, between the present code prediction and the numerical computation by 
Mei and Plotkin [11]. The agreement was within 2%. 

4. Also, the present code was used to compute an impulsively started flow over a backward-facing 
step with Re = 400 and a 1 : 2 expansion ratio. The solution was marched in time and the friction 
factor and reattachment length were compared with the steady flow results for the same case. 
The agreement was within 1%. Similar agreement was found upon examining a forward-facing 
step case. 

5. Finally, unsteady flow calculations were conducted for oscillating flow (zero mean) in a straight 
channel and were in excellent agreement with available analytical solutions for fully developed 
channel flow [2]. 


RESULTS AND DISCUSSION 

Table 1 lists the cases studied in this paper. These cases cover wide ranges of Re^ (from 187.5 
to 2000), Va (from 1 to 10.66), the expansion ratio (1:2 and 1 :4) and A r (2 and 4). Also, shown 
in the table are the three geometric cases discussed above: 
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Fig. 2. Time-varying sinusoidal velocity at the channel inlet. 


(a) Asymmetric expansion/contraction, cases A, D, E, J, K and L. 

(b) Symmetric expansion/contraction, cases B, F and H. 

(c) Symmetric blunt body, cases C, G and I. 

All cases correspond to operating conditions of the NASA SPRE heater at St = 0.021333. The inlet 
velocity at either end of the channel varies sinusoidally with time, as shown in Fig. 2. 

Figure 3 shows the envelope in which different Stirling engines operate, plotted in terms of Re^ 
vs Va [4]. In the figure, different criteria [14, 15] for the transition from laminar to turbulent flow 
are shown, for a straight channel. Below these lines (low Re^) the flow will remain laminar 
throughout the cycle, while above them some combination of laminar/transitional/turbulent flow 
occurs over the cycle. As evident from the plot, most of the Stirling engine conditions are in the 
transition or “fully turbulent” zone. Efforts are underway to map the conditions under which 
quasi-steady turbulence models can be applied to oscillating flow conditions. For more details, see 
Ref. [16]. 

The solid circles shown in Fig. 3 are for cases A, E and K (see Table 1) as well as other runs 
made (not shown in this paper). We have attempted to identify the region in which the flow is 



Fig. 3. Envelope in which different Stirling engines operate, together with: (i) the criterion for tran- 
sition from laminar to turbulent flow in straight channels, Re^„ * 822 x (Va a5 ), Re^ = 770 x (Va 0 * 5 ) and 
Re^a, = 453.23 x (Va 0 * 57 ); (ii) the different cases studied in the present work. 
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laminar quasi-steady. Case A lies in that region, while cases E and K are in the non-quasi-steady 
region (to be discussed in more detail later). The shaded area in the figure indicates our estimation 
(outcome of the present analysis) that below it the flow is laminar quasi-steady, while above it the 
flow becomes laminar non-quasi-steady. 

Effect of the Geometry 

Figure 4 shows the streamlines of the oscillating flow at different velocity phase angles of 30, 
60, 90, 120, 150, 180, 210, 240, 270, 300, 330 and 360. The results are for Re max = 187.5, Va = 1 
and an asymmetric expansion/contraction with a ratio of 1 : 2 (case A). Figure 5 shows similar 
results for a symmetric expansion/contraction (case B); while Fig. 6 is for a symmetric blunt body 
(case C). From the plots it can be seen that in the sudden expansion process, the separation bubble 
behind the step grows gradually during flow acceleration and then shrinks gradually during flow 



Fig. 4. Streamlines of oscillating flow at different velocity Fig. 5. Streamlines of oscillating flow at different velocity 

phase angles, for Re^ = 187.5 and Va = I. (Case A, asym- phase angles, for Re^ « 187.5 and Va = 1. (Case B, sym- 
metric expansion/contraction with a ratio of 1:2.) metric expansion/contraction with a ratio of 1:2.) 
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deceleration. However, this bubble disappears completely during the flow reversal (sudden 
contraction). Also, as expected, the size of the separation bubble (at a given velocity phase angle) 
gets progressively bigger from A to B to C. These cases, as described above, show a quasi-steady 
behavior. 

As for the friction factor results (not shown in this paper for all cases); it was found that for 
case B the friction factor is of the same order of magnitude as the asymmetric expansion (case A). 
The difference lies in the developing zone beyond the reattachment point. Similar trends are 
observed for Re max = 1000 — cases E, F and G. 

Effect of Re nw and Va 

Figure 7 shows the streamlines of the oscillating flow at different velocity phase angles of 30, 
60, 90, 120, 150, 180, 210, 240, 270, 300, 330 and 360. The results are for Re max = 500, Va - 10.66 
and an asymmetric expansion/contraction with a ratio of 1 :2. Similar results are shown in Fig. 8 
for Re max = 1000, Va = 10.66 and an asymmetric expansion/contraction with a ratio of 1 : 2 (case 






Fig. 6. Streamlines of oscillating flow at different velocity Fig. 7. Streamlines of oscillating flow at different velocity 

phase angles, for Re^ = 187.5 and Va - 1. (Case C, sym* phase angles, for Re^ = 500 and Va = 10.66 (Asymmetric 

metric blunt body with a ratio of 1:2.) expansion/contraction with a ratio of 1:2.) 





Oscillating flow in channels with changing cross section 



E). Experimental and analytical results for a straight channel[4], also see Fig. 3, indicate that the 
flow can be assumed laminar for the two cases. Due to the lack of experimental evidence for the 
channel with a sudden change in cross section, the flow is considered laminar. It was observed while 
watching the animation videos of the above cases that as the fluid accelerates during the sudden 
expansion phase, the separation bubble region grows in magnitude and physical size. This effect 
is observed irrespective of the magnitude of the flow Re. However, during the flow deceleration, 
this separation bubble loses momentum rapidly and disappears for cases A-D (low Re) but 
continues to grow while losing momentum for cases E-L (relatively higher Re). This observation, 
can be seen also by comparing Figs 4, 7 and 8. 

For the cases presented, Va is also increased with the Re to maintain a constant St corresponding 
to the actual Stirling engine. 






Fig. 8. Streamlines of oscillating flow at different velocity 
phase angles, for Re mas = 1000 and Va = 5.33 (Case E, 
asymmetric expansion/contraction with a ratio of 1:2.) 


Fig. 9. Streamlines of oscillating flow at different velocity 
phase angles, for Re^ = 187.5 and Va = 1. (Case D, asym- 
metric expansion/contraction with a ratio of 1:4.) 
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Fig. 10. Reattachment length vs time for oscillating flow at Re^^ 187.5 and Va=l. (Case A, 
asymmetric expansion/contraction with a ratio of 1:2.) 


It was found from the analysis in this work that as Va increases, the flow becomes non- 
quasi-steady at a lower Re^. This is shown in Fig. 3 by the shaded area, below which the flow 
can be considered quasi-steady. Also, this could lead to the conclusion that as Va increases the 
flow is likely to become turbulent at a much lower Re max than normally observed for straight 
channels. 


Effect of the Expansion Ratio 

Figure 9 shows the streamlines of the oscillating flow at different velocity phase angles of 30, 
60, 90, 120, 150, 180, 210, 240, 270, 300, 330 and 360. The results are for Re^ = 187.5, Va= 1 
and an asymmetric expansion/contraction with a ratio of 1:4 (case D). 

Comparing cases A and D (both are similar except for the expansion ratio), the recirculation 
zone in case D grows to almost 4 times the size in case A. The recirculation zone, being larger, 
dissipates its energy during the fluid deceleration and experiences some growth. This affects the 
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Fig. 11. Reattachment length vs time for oscillating flow at Re^ = 1000 and Va = 5.33. (Case E, 
asymmetric expansion/contraction with a ratio of 1 : 2.) 
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friction factor by more than an order of magnitude (not shown in this paper). Also, because of 
a larger recirculation zone, which results in good mixing, wall heat flux is reduced (not shown in 
this paper). Again, by comparing cases A and D, it can be seen that although case A is quasi-steady, 
case D is not because of the higher expansion ratio. 

For a higher Re, cases E and J (Re^ = 1000) and cases K and L (Re m „ = 2000), a similar trend 
regarding the size of the recirculation zone and the corresponding effect on friction and heat 
transfer is observed. 


Reattachment Length 

Figure 10 shows the reattachment length vs velocity phase angle for Re^ = 187.5, Va = 1 and 
an asymmetric expansion/contraction with a ratio of 1 : 2 (case A). Also, shown is the reattachment 
length, given by Armaly et al. [5], for steady flow at the instantaneous Re. The results of the present 
analysis show good agreement, indicating that the flow is quasi-steady as explained earlier. On the 
other hand, Fig. 1 1 shows a similar plot of the reattachment length vs velocity phase angle at a 
higher Ren^, = 1000, Va = 5.33 (case E). It can be seen from Fig. 1 1 that the reattachment length 
will grow initially at a lower rate than in the corresponding quasi-steady case, thereafter it will 


(a) 



(b) 



x/s 

Fig. 12. (a) Friction factor vs time for oscillating flow at Re^, = 187.5 and Va = 1. (Case A, asymmetric 
expansion/contraction with a ratio of 1 :2.). (b) Normalized friction factor vs time for oscillating flow at 
Re„„ = 187.5 and Va = 1. (Case A, asymmetric expansion/contraction with a ratio of 1:2.) 
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continue to grow at a higher rate and then disappear during flow reversal, i.e. case E is 
non-quasi-steady. 


Friction Factor 

Figure 12(a) shows the instantaneous friction factor vs dimensionless axial distance at different 
velocity phase angles, for Re m „ = 187.5, Va = 1 and an asymmetric expansion/contraction with a 
ratio of 1:2 (case A). A similar plot is shown in Fig. 12(b) but with the coordinate being the 
instantaneous friction factor divided by the steady-state friction factor at the instantaneous Re. The 
quasi-steady behavior is observed by having values of ///„ close to 1.0 in Fig. 12(b) at all times 
and most of the channel axial locations. 

Figure 13(a) shows the instantaneous friction factor vs dimensionless axial distance at differ- 
ent velocity phase angles, for Re^ = 1000, Va = 5.33 and an asymmetric expansion/contraction 
with a ratio of 1 :2 (case E). A similar plot is shown in Fig. 13(b) but with the coordinate being 
the instantaneous friction factor divided by the steady-state friction factor at the instantaneous 
Re. It can be seen from Fig. 13(b) that f/f a departs considerably from 1.0, particularly after the 
step, indicating a non-quasi-steady behavior. The friction factor can be a factor of 2 higher or 
lower than the steady flow values for X/S east from the reattachment location. This indicates 



XI S 


(b) 



Fig. 13. (a) Friction factor vs time for oscillating flow at Re^, = 1000 and Va = 5.33. (Case E, asymmetric 
expansion/contraction with a ratio of 1 : 2.). (b) Normalized friction factor vs time for oscillating flow at 
Re — = 1000 and Va = 5.33. (Case E, asymmetric expansion/contraction with a ratio of 1:2.) 
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that the steady-state friction factor correlations cannot be used in such applications. Work is 
underway to develop new correlations for the friction factor and pressure drops for oscillating flow 
conditions. 


CONCLUDING REMARKS 

In this paper the oscillating flow (zero mean) between two parallel plates with a sudden change 
in cross section was studied. The flow was assumed to be laminar incompressible with the inflow 
velocity uniform over the channel cross section but varying sinusoidally with time. Under these 
conditions, the fluid undergoes a sudden expansion in one-half of the cycle and a sudden 
contraction in the other half. The flow reversal, under such conditions, is caused not only by the 
flow oscillations but also by the sudden change in cross section. 

A computer code, CAST, developed by Peric and Scheuerer [12] has been modified to handle 
time-varying boundary conditions. The CAST code solves Navier-Stokes equations in 2-D using 
a finite volume method. The code has been validated by comparing its predictions with available 
computational, experimental and analytical data for straight channels and those with a sudden 
change in cross section. Good agreements were found for the friction factor, reattachment length 
and minimum and maximum velocities at different axial channel locations. 

The computations were extended to oscillating flow conditions. The cases examined emulate the 
operating parameters of the NASA SPRE (Space Power Research Engine). The Ren,, in cases 
presented here is chosen to be sufficiently low that the assumption of laminar flow holds true during 
the entire flow cycle. The cases examined are summarized in Table 1. 

In all cases examined, a separation zone appears during the sudden expansion and grows as the 
flow accelerates. This growth, however, depends on the Re mai . For low Re mu (= 187.5), the growth 
follows a quasi-steady behavior, while for the higher Re ra „ ( = 1 000) the growth is very rapid and 
the separation bubble continues growing during flow deceleration (non-quasi-steady behavior). 

A shaded area has been identified (see Fig. 3), below which the flow is laminar quasi-steady. This 
indicated that as Va increases the flow becomes non-quasi-steady at a lower Re^ . 

When the flow reverses, the fluid goes through a sudden contraction and the recirculation bubble 
from the previous half-cycle is swept back into the smaller section of the channel. 

Such flow behavior causes the instantaneous friction factor to deviate substantially from the 
steady-state friction factor for the same flow parameters. The friction factor can be a factor of 2 
higher or lower than the steady flow values for X/S east from the reattachment location. Work 
is underway to develop new correlations for the friction factor and pressure drops under oscillating 
flow conditions. 

Upon examining the effect of the channel expansion ratio, it was found that for the same Re 
and Va, increasing the expansion ratio increases the friction losses. This is consistent with the 
observations made by several researchers studying steady flows with a change in the channel cross 
section. 
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CHAPTER I 


INTRODUCTION 


Inside the heat exchangers of an operating Stirling engine are 
hundreds of tubes where heat is transferred with the working fluid under 
oscillating flow conditions. Presently, NASA is using one- dimensional 
performance codes to model and simulate this unsteady, oscillating flow 
towards the development of the free-piston Stirling engines. These 
performance codes utilize incompressible steady-flow heat transfer and 
friction factor correlations, while excluding the effects of oscillation and 
the compressibility of the working gas. Thus, the validity of using the 
steady-state correlations for the oscillating conditions is questionable 
[Tew,1987; Ibrahim et.al.,1990]. Confidence is also lacking in the predicted 
performance by these one-dimensional analyses. 

The present study investigates the complex oscillating flow conditions 
within the Stirling engine. The objectives of the current research are as 


follows: 


1) to develop a two-dimensional modeling code for compressible 


unsteady, periodic internal flows; 

2) to better understand the fluid flow and heat transfer 
phenomena under the pulsating/oscillating conditions and to 
compare the results with the steady flow conditions; 

3) to investigate the compressibility effect on the friction and the 
heat transfer; and 

4) to visualize the results on video with computer animation. 

The current computer code evolved from the codes developed earlier 
by Chilukuri [1979], Madavan [1981] and Chiu [1984], whom have used it to 
estimate incompressible steady flow. However, the current effort is to 
modify the computer code to accept unsteady compressible flows. 

Two types of unsteady flows will be examined in this thesis, and to 
avoid possible confusion, the following definitions will clarify the 
differences between the oscillating and pulsating flows: 

1) Oscillating - the flow fluctuates back and forth in the channel 
and the time-averaged velocity of this flow type is equal to 
zero. 

2) Pulsating - the flow is fluctuating around a steady driven 
flow and the time-averaged velocity of this flow type is non- 
zero. An example of the pulsating flow is the circulation of 


z 2 


blood in one's veins. 


Review of Past Research 


Steady Flow (Incompressible) : 

Numerous studies exist in the literature on the problem of steady 
incompressible flow through parallel plates with various boundary 
conditions, and most of the pioneer works were based on the boundary 
layer theory. Some of the researchers were Schlichting [I960], Siegel and 
Sparrow [1959], and Bhatti and Savery [1977]. 

In his "Boundary Layer Theory" book, Schlichting [1960] analytically 
resolved the boundary layer equations for flow between parallel plates, and 
his analysis was based on the Blasius solution technique for external flat 
plate flow. For heat transfer analysis, Siegel and Sparrow [1959] obtained 
a closed form approximated Nusselt number equation for simultaneously 
developing flow under the condition of uniform wall temperatures, while 
Bhatti and Savery [1977] solved the same problem except the heating they 
used was provided by the constant wall heat flux. 

Narang and Krishnamoorthy [1976] presented the pressure and 
velocity solutions at the entrance flow regions of the parallel plates for the 
low Reynolds number. They solved the Navier-Stokes equations by 
linearizing the inertia terms and observed overshoots in the developing 
velocity profiles which was not obtained by the boundary layer equation. 
Their results were compared with data from Schlichting [1960] and they 
concluded the boundary layer assumptions for the entrance region were 
only valid for the higher Reynolds number (Re > 200). They also observed 


the v-velocity at the entrance was as large as 30% of the inlet u-velocity 
for the low Reynolds number (Re = 2), where the v-velocity is assumed 
zero in the boundary layer theory. 

Some other researchers who solved the full Navier-Stokes equations 
were Brandt and Gilles [1966], Morihara and Cheng [1973], and Chen [1973]. 
Brandt and Gilles [1966] utilized the finite-difference approximation in their 
study of incompressible flow in a straight channel with the presence of a 
transverse magnetic field. Morihara and Cheng [1973] solved the Navier- 
Stokes equation by first eliminating the pressure terms followed by 
quasilinearizing the nonlinear terms in the resultant momentum equation. 
Finite-difference formulas combined with the Gaussian elimination solver 
was utilized to reach the solution. The Reynolds number that was 

considered, ranged from zero to 4000. 

In the entrance flow study by Chen [1973], he performed the 
momentum integral method to solve the Navier-Stokes equation and obtained 
approximated closed form solutions for the center-line velocity development 
and the average axial pressure drop. His results were compared and were 
in excellent agreement with the finite-element solution and the experimental 
study at the low Reynolds numbers for flows in circular tubes by Atkinson, 
et.al. [1969]. A comparison was also conducted with the finite- difference 
analysis by Friedmamn, et.al. [1968]. Velocity overshoot in the developing 
profiles were observed by all the investigators who applied the full Navier- 
Stokes equations. 

While the majority of the experimental data for steady internal flows 


* H- 


were performed with the circular tube, a few studies on the parallel plates 
were found. Heaton, et.al. [1964] performed an experimental study of 
simultaneously developing flow with constant uniform wall heat flux in 
Annulus which simulated the parallel plates channel. In a more recent 
study, Lou and Barton [1973] assembled a recirculation loop experimental 
apparatus with a rectangular channel test section that measured 7 inches 
wide by 1/4 inch deep and was 2 to 6 inches long. The fluids examined 
were water (Pr=5.7), 20% ethylene-glycol (Pr=10) and 72% ethylene-glycol 
(Pr=50). Their heat transfer results for constant wall temperature 
compared well with the theoretical data from Lou [1971]. 

Finally, an experimental investigation of the developing thermal 
entrance region with heating at the lower plate and cooling at the upper 
plate was done by Kamotani and Ostrach [1976]. 

Variable Property Flows ; 

Two types of variable property flows have been considered in the 
literature: thermally expandable (temperature dependent) and compressible 
(temperature and pressure dependent). An initial thermally expandable 
flow investigation was studied by Van Driest [1952]. Van Driest modeled 
the steady laminar, external air flow over a flat plate using the boundary 
layer method. He assumed the Prandtl number and the specific heat as 
constants while employing the Sutherland law to describe the fluid 
viscosity. Van Driest presented predictions for skin-friction and for the 
heat transfer coefficient as functions of the Reynolds number, the Mach 
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number and the wall to free-stream temperature ratio. 


For the thermally expandable internal flows, Bankston and McEligot 
[1969] utilized the superposition method on their study of the coolant 
channels of a nuclear reactor. The objective was to determine the wall 
temperature from an imposed constant wall heat flux by the variable 
property model and the incompressible model. They found that the 
variation of the fluid properties tends to shorten the thermal entry region 
when heating the fluid. 

Later, Bankston and McEligot [1970] extended their study by using 
the numerical method of finite control volume to solve the two-dimensional, 
coupled boundary layer equations. The study was done on the circular 
tube with constant heat flux; laminar and turbulent flows were considered. 
They showed that the variable property slightly affected the Nusselt 
number but the friction parameters were tripled when compared with 
incompressible results under heating conditions. The same analyses were 
performed for the parallel plates channel with constant wall temperature 
in Schade and McEligot [1971] paper. The study was performed in the 
wall/inlet temperature ratios that ranged from high heating (T w /T fl =10) to 
extreme cooling (T w /T 0 =O.l). Similar conclusions as the circular tube study 
were found; the effects of gas property variation on the laminar flow of 
air are slight to moderate for heat transfer and more severe for the 
friction. 

Herwig et.al. [1990] introduced an efficient approximation method to 
resolve the variable viscosity pipe flow caused by temperature changes. 
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The basic idea is to expand the viscosity in a Taylor series so that the 


Navier-Stokes and energy equations can be decoupled. By doing so, no 
additional nonlinearities are imposed by the temperature-dependent 
viscosity, hence, the set of equations that has to be solved is much simpler 
than the full set of coupled equations. 

In the case of compressible flows, most literature found dealt with 
the study of external aerodynamic design of planes or high speed flows 
over a flat plate. MacCormack [1976] examined the time-dependent, two- 
dimensional compressible Navier-Stokes equations at Mach number of two. 
He presented a method to solve for the interaction of a shock wave with 
the boundary layer on a flat plate. The method time-splits the governing 
equations into a hyperbolic part and a parabolic part. The hyperbolic part 
was solved by the explicit method while the parabolic part was solved by 
the implicit scheme. 

Other methods to resolve the high-speed compressible flow over a 
flat plate were introduced. Beam and Warming [1978] and Kwon et.al. [1988] 
applied the implicit finite- difference algorithm with a second-order- time 
accurate scheme which required data storage of two time levels. 
MacCormack [1982] introduced a new efficient unconditionally stable method 
that consisted of two stages to solution. The first stage used an explicit 
predictor-corrector finite- difference method while the second stage 
transformed the governing equations into an implicit form. Han [1983] then 
modified the SIMPLE (Semi-Implicit-Method for Pressure-Linked-Equation) 
Algorithm to predict transient analysis of both high speed (compressible) 
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and low speed (incompressible) fluid flows. 

Besides flows over a flat plate, Kim [1990] studied the transonic flow 
over an axisymmetric curved hill. Hassankhan [1983] examined the effects 
of fluid compressibility on the stagnation point of a circular cylinder in 
cross flow. Hassankhan found that the friction increases when the Mach 
number and/or the wall/inlet temperature ratio increases. 

As for internal compressible flow, no literature was found for the 
parallel plates geometry but two studies were located for circular tubes. 
Presler [1971] studied helium gas flow in a circular tube with uniform 
entrance conditions and wall heat flux conditions. The study was done 
both analytically and experimentally with an electrically heated small- 
diameter Inconel tube. The numerical results were performed for low 
values of the initial Mach number to prevent the choking effects. Presler 
reported that for any level of uniform wall heat flux, the local Nusselt 
number showed only small deviations from the constant property analyses, 
while the friction was larger for higher wall heat fluxes. The comparison 
of his analytical and experimental results were all within 10 percent. 
Similar results were found by Cebeci and Bradshaw [1984]. 


CHAPTER n 


ANALYSIS 


Figure 2.1 reveals the cross sectional view of the NASA Stirling 
Power Research Engine (SPRE). The components of interest are the two 
heat exchangers (heater and cooler) that supply and eject the necessary 
thermal energy to continue the oscillating motion of the working helium 
gas. A closer look at the inside of these heat exchangers reveals 
hundreds of straight, narrow tubes that contain the oscillating helium gas. 
Instead of examining the global heat exchanger problems, the present study 
is to explore and understand the detail flow and heat transfer phenomena 
of the fluid within each tube. The parallel plate channel is chosen for 
modeling this detailed single tube problem. 





Geometry and Coordinate System 


Figure 2.2 shows the parallel plate channel geometry and the 
Cartesian coordinate system used. The origin of the coordinate system 
resides at the west corner of the lower plate with the x-axis set along the 
stream-wise direction while the y-axis starts from zero at the lower plate 
to the value of H/2 at the channel's center-line. For the steady and 
pulsating flows, the fluid enters the channel from the west where the 
numerical calculation begins (Fig. 2.3). As for the oscillating flow, the 
fluid enters from the west during the forward half of the cycle and then 
shifts to enter from the east during the reverse half of the cycle (Fig. 
2.4). 


Assumptions 

In the heater and cooler tubes, the unsteady periodic fluid flow 
contains the complexities of 3-dimensional, laminar, transient and turbulent 
flows; the fluid enters the heater and cooler tubes from the expansion and 
compression space, respectively, through the passage of a sudden change 
in the cross-sectional area; the oscillating phase angle difference between 
the displacer and the power piston causes real gas compression and 
expansion within the tubes. Due to the enormous amount of experimental 
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Figure 2. A: The boundary conditions for oscillating flows. 
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data and time it would take to model these complex conditions, the present 
study is simplified and is limited to the following assumptions: 

a) The fluid flow is a two-dimensional, laminar with variable 
properties except the and Pr are assumed constant; 

b) The inlet velocity profile is uniform and 

i. constant with time for steady flow; 

ii. varies sinusoidally with time for unsteady flows; 

c) The fluid enters the channel with uniform temperature and 
heat is transferred from (or to) the constant temperature plate 
walls; 

d) The axial viscous diffusion and heat conduction are negligible; 

e) The flow is subsonic. 


Static Temperature versus Stagnation Temperature 

To prevent confusion when examining compressible flow, the 
difference between static and stagnation temperatures will be clarified. 
The stagnation enthalpy of the fluid represents the combination of thermal 
energy and kinetic energy. Since the specific heat of the fluid is assumed 
to be constant, the stagnation enthalpy equation can be written in terms 
of the stagnation temperature (sometimes called total temperature): 


* /3 


,, j + „2 where: x is the stagnation teoperatuze , _ , . 

x b T + . - , , (2.1) 

2 <- p T 1b the static temperature 

As suggested by the equation, the stagnation temperature represents 
the total temperature as the gas is brought to rest by an adiabatic 
process. If the fluid velocity is slow or the Mach number is close to zero, 
the fluid kinetic energy is negligibly small compared to the fluid thermal 
energy; hence, the stagnation temperature and static temperature are 
equal. 

H » 0 - T ■ T 

For the incompressible and thermally expandable flows, the 
presumption of the zero Mach number results in identical static and 
stagnation temperature values, therefore, no distinction between these 
temperatures is necessary as in most incompressible and thermally 
expandable flow literature. However, the compressible flow presumes a 
non-zero value of the Mach number which means the kinetic energy of the 
fluid can no longer be ignored and the two temperatures are differentiated 
by equation 2.1. 
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Governing Equations 


The complete Navier- Stokes (NS) equations can model any viscous 
flow problems. Mathematically, NS equations are elliptic partial-differential 
equations which in a physical view-point, the flow at any location is 
affected by the rest of the flow field. Computing these coupled elliptic 
differential equations by finite- difference methods requires iterative 
solvers, but these solvers require heavy computer time and memory space. 

For flows with a predominant flow direction, as in the current study, 
where the axial viscous diffusion and heat conduction are small and 
negligible, the NS equations can be simplified into the Partially-Parabolized 
Navier-Stokes (PPNS) equations. The PPNS equations are parabolic partial^ 
differential equations by nature except the elliptic behavior associated with 
the pressure remains. Having the streamwise diffusion neglected, the flow 
at any location is affected only by the upstream flow and can be solved 
by marching methods which are relatively more efficient in computer time 
and memory space than iterative methods. Please note that only the 
pressure in the PPNS equations is computed by using iterative solvers. 

Utilizing the assumptions stated earlier, the governing PPNS 
equations in primitive variables are as follows: 


Continuity Equation: 



X-Momentum Equation: 



Y-Momentum Equation: 



Note that the viscous diffusion energy is accounted for in the last two 
terms of the energy equation, hence the static temperature is used as a 
primitive variable. 

The pressure has an elliptical effect throughout the flow field and 
is described by the Poisson equation [Anderson et.al,1984] which can be 
derived from the momentum equations. By rearranging equations 2.3 & 2.4, 
the pressure gradients can be written in the form: 



( 2 . 6 ) 






(2.7) 
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and the Poisson equation is: 


V,p *B*f£’TF*Tp (2 ' 8) 

The above governing equations are similar to the equations Chiu 
[1984] utilized in his analyses, except Chiu examined the steady, 
incompressible flows which made the apparent form of his PPNS equations 
somewhat further simplified. The present equations are more general as 
they contain additional time dependent terms to account for the unsteady 
flows and the equations have variable fluid properties (density, viscosity 
and conductivity) within the derivatives. 


Equations of State 

To account for fluid compressibility, one state equation is required 
for describing each variable fluid property (density, viscosity and 
conductivity). Assuming the gas is ideal, the density can be determined 
from the pressure and the static temperature by using the Perfect Gas 
Law. 

p - -£=. where R la the Gas constant (2.9) 

r RT 
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To date, the viscosity of gases has been described by two commonly- 


used equations, the Power Law and the Sutherland Law. 


The Viscosity Power Law: 



The Sutherland Law: 


where a is a constant 


( 2 . 10 ) 



where B and S are Sutherland coefficients 


( 2 . 11 ) 


The a, B and S are constants depending upon the gas used (their values 
for air and helium are listed in Table 2.1). The Power Law and the 
Sutherland Law relate the viscosity with the temperature of the gas. From 
the experimental data by Tribus et.al. [1942] for viscosity of air at various 
temperature. Driest [1952] found that the predictions from the Sutherland 
Law matched the data more closely than the predictions from the Power 
Law. Hence, the Sutherland relationship was chosen for this study to 
achieve more accurate viscosity estimations. 

As for the fluid conductivity, again, two relationships were commonly 
used; the conductivity Power Law and the definition of the Prandtl number. 


The Conductivity Power Law: 



where 0 is a constant 


( 2 . 12 ) 
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The Definition of the Prandtl Number: 


k - 



(2.13) 


For most gases, the Prandtl number and specific heat are relatively 
constant for a wide range of temperature which makes the gas conductivity 
directly proportional to the gas viscosity when incorporating the Prandtl 
number definition. In this study, the Prandtl number definition is applied. 
As these state equations are normalized, it will become apparent why the 
Prandtl number is more advantageous over the Power law. 


Boundary Conditions 

As shown on Figure 2.2, the computational domain is one-half of the 
channel, from the bottom plate to the center-line of the channel. At the 
lower solid wall boundary, the no-slip wall condition applies and the wall 
temperature is constant. 

y = 0 — u « v » 0 . T = T v 


20 


While at the line of symmetry, u-velocity and temperature gradients are 
assumed to be zero. 


y - 


B 

2 


v - 0 , 


du 

■ay 


0 , 


dT 

Sy 


- 0 


The upper and lower boundary conditions are identical for the steady, 
pulsating and oscillating flows. 

Ensuring the incoming fluid contains the same energy level, the 
imposed inlet stagnation temperature remains constant for all flow types, 
while the inlet v-velocity is assumed to be zero. 

t 0 * constant , v„ * 0 


The imposed inlet u-velocity, however, differs according to the flow 
type. For the steady flow, inlet u-velocity is constant (see figure 2.3). 

u„( t) = u„ = constant 


In the case of pulsating flow, the u-velocity varies sinusoidally around a 
mean flow according to the relation: 

U e {t) “ Uw*mn.o' (1 + % Sin<*t) 

(2.14) 

where: is the imposed velocity fluctuation amplitude 
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Finally, the imposed inlet u- velocity for oscillating flow varies around 
a zero-velocity. 

u c (t) « Sin at (2.15) 

Referring to figure 2.4, the inlet flow enters the channel from the west 
when u fl is positive, and east when u fl is negative. 

The boundary conditions for the pressure Poisson equation will be 
discussed in the next chapter. 


Similarity Parameters 


The following characteristic parameters are used to generally 
describe the fluid flow and heat transfer in similar systems. For 
incompressible and steady flow, only the Reynolds number and Prandtl 
number is needed for generalizing similar systems. 

Reynolds number (non-dimensional mass flux): 





(2.16) 


Prandtl number (ratio of viscosity over conductivity of the fluid): 


Pz « 


*x«f 


(2.17) 
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The additional Mach number is required when dealing with 
compressible flow. 

Mach number (non-dimensional flow speed): 

M zat ~ where : C is the speed of sound ( 2 . 18 ) 

C zof 


For unsteady flow, the Valensi number generalizes the fluctuating 
frequency of the system. 


Valensi number (non-dimensional frequency): 


Ci> 
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( 2 . 19 ) 


An amplitude ratio is needed to describe the geometric similarity 


between oscillating flow systems. 


Amplitude ratio (non-dimensional fluid displacement): 


* ~ 2 L 

( 2 . 20 ) 

where: A R < 1 -* fluid oscillates without exiting the channel 

Ag > 1 ■* fluid traverses quickly via the channel 
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Nondimensionalization 


The advantage for nondimensionalizing the variables is that the 
characteristic parameters such as the Reynolds number, the Prandtl 
number, the Mach number, etc., can be varied independently. Also the 
values of the normalized flow variables often fall in the range such as zero 
to one. 

Location of the reference values depends on the flow type. The 
subscript "ref" used throughout the thesis has the following definitions: 

a) For steady flow: ref -* the value at the channel's inlet. 

b) For pulsating flow: ref -» the mean value at the channel's 

inlet. 

c) For oscillating flow: ref ■» the maximum value at the 

channel's inlet. 

With the appropriate reference values, all the variables are 
nondimensionalized as follows: 


* 


£T = — H- 


the normalized streamrise velocity 


the normalized transverse velocity 


Pr.f u z,l 


the normalized pressure 


_ T 


the normalized static temperature 


the normalized stagnation temperature 


the normalized x distance 


Prof Ur.fy 


the normalized y distance 


Hr./ 


the normalized time 


the normalized fluid density 


the normalized dynamic viscosity 


.. . * 


the normalized fluid conductivity 


the normalized streamvise mass flux 


V* - p* V 


the normalized transverse mass flux 
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The Normalized Governing Equations 


Thus, the partially-parabolized Navier-Stokes equations (2.2— 2.5) can 
be converted as follows: 

Continuity Equation: 



( 2 . 21 ) 


X-Momentum Equation: 


dU dU dU m „dP + 3 dU\ 

p ~3x "Syr Sy) 


Y-Momentum Equation: 


,'#..rg.r£ 


_3P + 4 3 /„. dv\ 
ly 3 ‘3?r 75?) 


Energy Equation: 


• •98 . rr» 36 . » t 9 

p Sx*^ 


38 _ 1 3 /l.. 38\ 
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( 2 . 22 ) 


(2.23) 


For the Poisson equation, the x and y-momentura equations are rearranged. 

£ ■ +•£*«•£*»•£-*(*•£)] - » «*•»> 

< 2 - 2 ‘> 


The pressure Poisson Equation: 


Note that for thermally expandable flows, hence the kinetic 

energy terms are negligible and the energy equation can be simplified to 
balance only the thermal energy. As for incompressible flows, p\ u\ k* = 
1 which also converts U* = U and V* = V, and the coefficient value 4/3 is 
changed to one in the y-momentum diffusion term due to simplification by 
the incompressible continuity equation. With these reductions and for the 
steady flow condition, the normalized governing equations will be identical 
to the equations Chiu [1984] utilized. The present governing equations are 
more general than Chiu’s, and they are capable to describe the different 
flow types considered in the present analyses. 


Normalized Equations of State 


Density (Ideal Gas Law): 


P* - 



(2.28) 


vheret a P ■ the change of P from the channel 'a inlet. 

For subsonic flow where « 1, the second term on the right hand 
side of equation 2.28 will be insignificant compared to the first term; 
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hence, for low Mach number flow, the second term can be ignored and the 
value of the density depends solely by the temperature only. This type 
of flow is assumed to be thermally expandable. 



The definition of Prandtl number is used because it equates the 
normalized viscosity and the normalized conductivity, and there is one less 
variable to solve. 
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Normalized Boundary Conditions 


The nondimensional boundary conditions are 
a) U * V » 0 , 0 « 8„ no-allp boundary conditions at the wall 


b) V - 0 


du 

■ 3 ? 


00 

3? 


symmetry boundary conditions 


c) T 0 » constant , V 0 » 0 imposed inlet temperature and v-velocity 

d) imposed inlet u-velocity 


i) 

u 0 {*) 

- a o -l 

for steady flow 

(2.31) 

ii) 

a o (*) 

■‘•HO 

for pulsating flow 

(2.32) 

Hi) 

U 0 C) 


for oscillating flow 

(2.33) 


Again, the normalized pressure boundary conditions will be discussed in 
Chapter III. 
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CHAPTER m 
METHOD OF SOLUTION 


Based upon the method of solution used in Chiu's thesis [1984], the 
present solution techniques have been revised and extended to cover a 
wider range of fluid flow problems including the ability to resolve 
unsteady periodic flows and compressible fluid flows. This chapter 
consists of six sections. The first section describes the staggered, finite- 
difference grid u tiliz ed in the computation domain. Sections 2 to 5 present 
the various finite-difference schemes applied on the governing equations 
and the numerical solvers used. Section 6 summarizes the chapter by 
illustrating the flow chart and the steps proceeded by the present 
computer program to reach the results. 
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The Finite-Difference Grid 


Having selected the Cartesian coordinate system for the computation 
domain, the present finite difference analysis uses a staggered grid 

I 

[Patankar,1980] to prevent the probable development of a wavelike pressure j 

or a wavelike velocity field. The idea behind the staggered grid is to 
define three different grids; two for the velocity components and one for 
the scaler variables. As shown in Figure 3.1, the solid dots indicate the 
grid locations for the scaler variables (the fluid properties, pressure, 
temperature and the velocity correction potential 4>), whereas the velocity 
components, indicated by the arrows, reside midway between these dot 
points. The horizontal arrows denote the locations for U and the axial 
pressure gradient FI, while the vertical arrows designate the positions for 
V and the transverse pressure gradient F2. Figure 3.1 also reveals the 
spatial step variables used on the unequally spaced grid. 

Despite the fact that the variables are actually defined at different 
locations, a single reference set of grid indices is used for convenience. 

Referring back to Figure 3.1, the label (i+1, j) actually identifies the 
spatial locations of the two arrows and one solid point enclosed by the 
dashed boomerang, thus, j is to the right of j and j is 
beneath P i+ ^ j, etc. 

Unlike the convectional grid which only requires a single control 
volume at each node, the staggered grid requires three distinctive spatial 
control volumes (Fig. 3.2) for the different governing equations. This 
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Figure 3.2: The three distinctive spatial control volumes for: 

a) x-momentum equation, 

b) y-momentum equation, 

c) continuity, energy, and pressure Poisson equations. 



concept leads to the different locations of the expansion point according 
to the equation; care is called upon when differencing the governing 
equations to the corresponding control volumes. For example, the terms 
in the x-momentum equation are expanded at the (i+1, j) u-velocity arrow 
(Fig. 3.2a), while the terms in the energy equations are expanded at the 
(i+1, j) pressure point (Fig. 3.2c). 

Figure 3.3 shows the computational domain of the parallel plates 
combined with the staggered grid. Notice that a set of fictitious points is 
positioned below the solid wall (j = 1) and another set is located above the 
line of symmetry (j = NJ+1); these fictitious points are used in describing 
the boundaries of the domain because the horizontal boundaries in 
staggered grid are placed along the v-velocity arrows. Hence, the no-slip 
and constant temperature boundary conditions stated in Chapter Two 


become: 


Vi. a ■ -Vi.i 

(3.1) 

V l.» - 0 

(3.2) 

e * » m e 

2 

(3.3) 


At the upper boundary, the symmetric boundary conditions are prescribed 
as: 

Vl.ai * ^l.ar » i (3.4) 

Vi.ms+i * 0 (3.5) 
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locations of boundaries. 






(3.6) 
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= e 


i.« 7*1 


The velocity and temperature profiles at the inlet boundary are 
prescribed at the same grid indie of i=l, but the actual locations of the 
profiles are slightly different due to the nature of the staggered grid. 

Figure 3.4 presents the two types of grid set-ups used in the 
present study. For the steady flow and the pulsatile flow analyses, the 
axial grid spacing is fine at the channel's entrance then increases down 
the stream while the transverse steps are fine closer to the lower wall but 
coarse at the channel's center (Fig. 3.4a). This set-up provides sufficient 
node points at the critical regions, like areas close to the entrance and the 
wall, where rapid gradient changes occur. As for the oscillating flow 
analysis, the set-up grid shown in Figure 3.4b is used and is similar to 
the steady flow grid, except for the fine axial steps that appear at both 
ends of the channel. The reasoning behind this is that during half the 
cycle, the flow is moving forward while the second half of the flow is 
moving backwards. 


Finite-Difference Formulation of the Momentum Equations 


Coefficients of the convective terms: 


The momentum equations are algebraically nonlinear differential 

4 Z6 




Figure 3.4: The x and y grid formats used for: 

a) steady flows and pulsating flows. 

b) oscillating flows. 
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equations. Solving these equations by numerical method requires 
linearizing the equations and computing them by an iterative procedure. 
The most common linearization strategy is done by lagging or extrapolating 
the coefficient U* and V* in the convective terms of the momentum 
equations. In this thesis, the coefficients in the x-momentum and y- 
momentum equations are denoted by the subscripts x and y, respectively. 


Lagging the coefficients: The coefficients are evaluated from the 

previous iteration values. This method is used where values of velocity 
components are stored. 
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Extrapolating the coefficients (regular): The coefficients are 

extrapolated by the upstream values. This method is always used in the 
first "Global” iteration at each time step, due to the lack of previously 


calculated values. 
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Extrapolating the coefficients (special case): For regions where only 

one upstream station is available (the second x-station in the channel's 
inlet), the coefficients are evaluated by the previous station. 
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Convective terms in the x-direction (axial) : 
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du 

Sx 


u m 


dV 

dx 


Except at the entrance regions with only one station in the upwind 


direction, the two convective terms are differenced by the three-point, 
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second-order accurate, upwind scheme. The polarity of the coefficient U* 
determines the direction of the upwind scheme (forward upwind for 
positive U , while backward upwind for negative U*). Also, the second 
order upwind differencing used in the present work was tested for the 
numerical diffusion and compared with a first order upwind differencing 
(normally used in the literature); the present scheme showed an 
improvement in the predictions of about 10%. 

For the x-momentum equation with forward going flow, the convective term 
becomes: 
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For the y-momentum equation with forward going flow, the convective term 
becomes: 


U ‘5x) i . U j I S* 


( 3 . 20 ) 


T _e*i,n*i 

i.j + Vi * ul ~ LLhi 
AX' AX * AX' 


4 40 


However, for the inflow region where only one station is available in 
the upwind direction, a two-point, first-order accurate, upwind scheme is 
used instead of the three-point scheme. 


For the x-momentum equation: 



1+UJ 


—0 t+Un+1 


t+l,n+l rn # a*l 

U 1+Uj ~ U l.j 

AX„ 


( 3 . 21 ) 


For the y-momentum equation: 

. „t*X,n»X 
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When the coefficient U* of the convective terms is negative, the flow 
is reversed and the forward moving flow differencing scheme is no longer 
valid. For this type of flow, a three-point upwind differencing scheme is 
still utilized except the differencing direction is reversed in the x- 
direction. 


Thus, for the x-momentum equation: 
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For the y-momentum equation: 
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(3.24) 


While the forward going differencing scheme uses velocities at the 
current n+1 iteration level, the reversed differencing scheme makes use of 
U u2 j * u i* 3 , j » U {*2 j and U i+ j j velocities at the previous n iteration level. 
The reason for using previous velocity values is because they have not 
been updated by the marching process at the current n+1 iteration level. 
Thus, storage for the n level velocities is required for the flow reversal 
region. 


Convective terms in the v-direction (transverse) : 


dY 


V 


dV 

dY 


The hybrid differencing scheme [Anderson et.al. 1984] is utilized to 
avoid numerical instability that a pure central differencing may resolve at 
a high mesh Reynolds number. A pure upwind differencing scheme may be 
used but the numerical diffusion may become too large compared to the 
actual diffusion at a low mesh Reynolds number. A scheme that switches 


♦ 


from central to upwind after a critical mesh Reynolds number may also be 
used but the sudden change could distort the convergence to a solution. 
To avoid this sudden switch, the present hybrid scheme applies a weighted 
average of the two schemes for large mesh Reynolds numbers and reduces 
to a pure central differencing scheme when the mesh Reynolds is smaller 
than the set critical mesh Reynolds value, R ( . This R c value should be less 
than or equal to 2.0 as suggested by Spalding [1972]. In the present 
analyses, R c is set to 1.9 (the same value used by Chiu [1984]) and no 
attempt was made to fine tune this number in maximizing the stability nor 
for minimizing the numerical false diffusion. 


For the x-momentum equation, the convective term becomes: 
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where W, A, and B are determined as follows: 
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R e * Critical mesh Reynolds number 
If Rl > Re . then A-l, B-0, and W-R^/R^ 
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If Ri i Re and & i -.Re * CAen A-0, B=0 , and AT= 1 


For the y-momenturo equation, the convective term becomes: 
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where the mesh Reynolds numbers are computed as follows: 


• t+l.B+l 

d* . y i*i,j . a v’ 


y yj+\.i 


D~ . i*l.i . A V* 
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By using these mesh Reynolds numbers, the W, A, and B variables can 
be determined by the same criteria as used for the x-momentum convective 
term. 


Pressure gradient anri diffusion terms ; 


dp 

Jx • 
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A second-order accurate, central differencing is used for these 
terms. Note that because the staggered grid is used, the expansion points 
for the central differencing equations are located at (i+1, j) where the 
velocities are evaluated (Fig. 3.2a, b). For the x-momentum equation, 

^ t+i,n _ t ♦x,n 

_ P - fj.i.j (3.27) 
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and for the y-raomentum equation, 
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For points immediately within the solid wall, the finite-difference 
representation of , equation 3.28 can be a poor representation 

of the diffusion term due to the use of a fictitious point outside of the 
boundary. A better representation is to assign a point at the boundary 


instead. For points immediately above the lower wall. 


^(ti'!?) becomes: 
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Time-dependent (unsteady! terms : 



P 


. dV 


A first order accurate in time, backward differencing scheme is 
chosen for the time dependent terms. Any finite- difference scheme would 
require storage of all values of parameters evaluated at earlier time steps. 
While a second order accurate scheme could be used, such a scheme would 
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require values of parameters at the two previous time steps to be stored 
which could require excessive amount of computer memory. Whereas using 
a first-order accurate scheme only requires the values of the previous time 
step, will result in saving half the required memory size. 


For the x-momentum: 
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For the y-momentum: 



v i* i.j 
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Momentum equations in Thomas Algorithm formats : 

In substituting the corresponding finite-differenced terms back into 
the momentum equations 2.22 and 2.23, the equations can be written for 
(i+1, j) grid point and can be rearranged to the following formats. 

The x-momentum: 
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for forward flow 
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The y-momentum: 
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for forward flow 
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The equations 3.34 and 3.42, when written for each y-grid point at . 
the i+1 axial station will result in a tridiagonal coefficient matrix for each 
equation. Combined with the proper boundary conditions (eq. 3.1, 3.2, 3.4 
& 3.5), the momentum equations are resolved by the efficient Thomas 
Algorithm solver. 

Points like the second stations and points adjacent to a solid 
boundary can be similarly constructed. 


Finite-Difference Formulation of the Continuity Equation 

If the exact pressure and temperature fields are used in the 
momentum equations, the resulting velocities will automatically satisfy the 
continuity equation. However the calculation procedure begins with 
estimated pressure and temperature, thus the velocity solutions from the 
momentum equations will not balance the continuity equation. An iterative 
procedure is necessary to reach the correct pressure and temperature 
fields. The present solution procedure adjusts the tentative velocities from 
the momentum equations by the continuity equation, and the corrected 
velocities are used to update the pressure (using the Poisson equation) 
and the temperature (using the energy equation). A detail description of 
the steady, incompressible velocity correction procedure can be found in 
Chiu [1984] and in Anderson et.al. [1984]; however, the present velocity 
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correction procedure is generalized to also work on compressible and 
unsteady flows. 

The continuity equation (eq. 2.21) is converted by the central 
differencing scheme except the unsteady, time dependent term is 
differenced by a first-order backward scheme. The expansion point is 
located at the (i+1, j) dot point. For the exact velocities values (U g * and 
V g *) and the exact density values (p e *), the continuity equation in finite- 
difference form becomes: 
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These exact values can be defined as the sum of the provisional values and 


the correction values. 
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. e*i.a*i . . e*l.a.l (3.52) 

v i'l.J * 


P*i*i,J = P i*l.J + Pci »! .3 


(3.53) 


The U*, V* and p* are provisional solutions obtained from the momentum 
equations, and the ideal gas equation, while U c *, V* and p c * are correction 
values. Substituting these definitions into equation 3.50 and rearranging 
all the correction variables to one side of the equation: 
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where is known since U , V and p have already been calculated. 

Equation 3.54 contains three unknown variables (U c *, V * and p c *) and 
in order to resolve this equation, further simplifications and variable 
reductions are necessary. 


This condition is true because the velocity at the previous axial 
location has already satisfied the continuity equation. 


b) Izzotational velocity cozzections 

This assumption permits the use of a velocity correction potential to 
relate U c * and V *, such that 
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this implies = 0 


c) 


. t+i.n+l 
<P i+2,j 


= 0 


This assumption is true as convergence is achieved. 


,, .c*l,n*l »e „ 

a) = 0 * “ 0 

These density correction terms will equal to zero when convergence 
is achieved. 


The above four assumptions simplify equation 3.54 and can be 
converted in the form: 
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As for the boundary conditions, the no-slip condition on the solid 
wall and the symmetric condition at the center-line are utilized. 
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Ar 
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where j = 2 or NJ+1 


The Thomas Algorithm solver is used to calculate this tridiagonal matrix. 
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Finite-Difference Formulation of the Pressure Poisson Equation 


Central differencing is used for the Poisson equation with the 
expansion point located at the (i+1, j) pressure point. 
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The pressure source term, Sp, can be determined by the FI and F2 
pressure gradients resulting from the same finite- difference forms of the 
momentum equations (eq. 2.25 and 2.26) using the corrected velocities from 
the continuity equation. Combined with the Neumann boundary conditions. 
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the resulting set of equations for pressure is solved by the method of 
Successive Over-Relaxation (SOR) by points. The use of an over-relaxation 
factor, usually in the range of 1.5 to 1.99, was necessary to efficiently 




reach the pressure solution with minimum pressure iterations. No in depth 
study was conducted on determining the best over-relaxation factor to use, 
but a trial and error process has shown that for unsteady flow and 
compressible flow, a minimum value of 1.9 is required for rapid 
convergence compared to the value of 1.5 used for most incompressible flow 
types. Moreover, the pressure source term (Sp) must be under- relaxed. 
Normally this under-relaxation factor starts with a small value (about .01) 
and gradually increases with each "Global" iteration (see the Solution 
Procedure Section) with a maximum value of .50 used. For variable 
properties flow, this under-relaxation factor generally remains below .30 for 
the best stable convergence. 


Finite-Difference Formulation of the Energy Equation 

In the present study, the momentum equations and the continuity 
equation are solved prior to the energy equation. The expansion point for 
the energy equation is located at the (i+1, j) temperature point. The terms 
in the energy equation are differenced as follows: 
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Three point upwind scheme 
Hybrid scheme 
Central differencing 
Backward differencing 


The followings are finite-difference expressions for the terms in the energy 
equation: 
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where W, A, and B are determined as follows: 
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Again, the value of R c is set equal to 1.9 for the present 
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By substituting equations 3.61 through 3.69 into equation 2.24, the 


energy equation can be rearranged in the form: 
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Combined with the constant wall temperature boundary condition (eq. 
3.3) and the upper symmetric boundary condition (eq. 3.6), the resulting 
Tridiagonal matrix is again solved by the efficient Thomas Algorithm 
method. However, for locations such as the 2nd x-station and points 
adjacent to the boundary walls, a 1st order upwind scheme for the axial 
convective term and a more accurate scheme for the thermal diffusion 
terms are used instead. 


Solution Procedure 


Based upon the Partially Parabolized Navier-Stokes (PPNS) equations 
solution procedure for incompressible, steady, subsonic flows [Anderson 


as zg, 


et.al.,1984], the present technique is extended to cover unsteady. 


compressible flow problems. For the steady flow calculation, the solution 
procedure found in Chiu's [1984] thesis is used with the energy equation 
resolved at each "Global" iteration (see Fig. 3.5) for variable fluid 
properties. As for unsteady flows, an extra "Time Step" loop is required 
to calculate the time dependent terms and to advance time increments in 
the fluctuation cycle. A flow chart of the present computer code is 
provided in Figure 3.5 and the computational steps are explained in the 
following summaries. 


Steady flow solution procedure summary ; 

51) Read all the input information, nondimensionalize all the 

variables, and initialize all counters. The normalized fluid 
properties are set temporarily to 1.0 while the initial boundary 
conditions at the channel's inlet are determined. 

52) The "Global" iteration level is advanced by one to n+1 level. 

53) Advance the axial counter "i" by one to the next stream-wise 
station. Use the n iteration level results or initial guessed 
values (for the first ’’Global" iteration) of the pressure and 
fluid properties, and determine the tentative U and V velocities 
by evaluating the momentum equations with the appropriate 
boundary conditions. 

54) Next, the continuity equation is balanced to correct the 
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Figure 3.5: Computational flow chart of the present solution procedure for 

the steady and unsteady fluid flows. 


to 


TIME STEP LOOP 















provisional U and V velocities. This requires the 
transformation of the correction velocities into the correction 
velocity potential 4> and solving equation 3.57. The residual 
mass source term is determined. 

55) Utilizing the corrected U and V velocities at the n+1 iteration 
level, the pressure source term Sp is evaluated at the i+1 axial 
station. 

56) For problems involving transfer of heat, the revision of the 
thermal field is performed. The temperature is updated to the 
n+1 iteration level via solving the energy equation with the 
latest values of U and V. 

57) Steps (S3) to (S6) are repeated until the end of the channel 
is reached. The U and V velocities, the temperature, and the 
Sp are revised for the entire flow channel to the n+1 "Global" 
iteration level. 

58) The pressure of the entire flow field is updated to n+1 
iteration level by utilizing the Poisson equation with the 
previously calculated Sp values and the Neumann boundary 
conditions. As mentioned earlier, this pressure correction 
process uses the Successive Over- Relaxation (SOR) by points 
method. The SOR method is an iterative procedure which 
sweeps the entire channel for the number of prescribed 
iterations. 
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S9) The last step in this n+1 "Global" iteration corrects the fluid 
properties by computing the equations of state with the latest 
temperature and pressure. 

S10) This "Global" iteration process repeats from step (S2) to (S9) 
until the residual mass ratio is converged to an acceptable, 
imposed percentage value (usually .1%) or until the number of 
"Global" iterations exceed the maximum prescribed value. In 
either case, the program prints out the detail flow results into 
the appropriate output tiles. A warning message will appear 
in these files if the program exit was due to an insufficient 
number of imposed iterations. 


Unsteady flow solution procedure summary : 

The unsteady flow solution procedure is based upon the steady flow 
solution steps with an additional time loop. At each time step, the 
unsteady PPNS equations are solved in a similar manner as in the steady 
flow procedure. The unsteady flow procedure is slightly different for 
pulsatile flow from the oscillating flow. For pulsatile flow, the fluid always 
travels in one direction, hence the numerical marching sweep moves from 
left to right. As for the oscillating flow, the fluid moves forward for half 
the fluctuating cycle and then it moves backward in the second half. 
During the backward flow period, instead of marching in the reversed 
direction, the present code "flips" the channel around as the input axial 
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velocity direction is changed. This actually simulates the same reversed 
marching but, this way, the code maintains the left to right sweeping 
direction throughout the oscillating cycle. The following outline briefly 
describes the unsteady flow solution procedure: 

Ul) The program reads all input information and sets up the "time 
step" loop. The program also determines if the flow is 
pulsatile or oscillating and sets up the appropriate inlet 
velocity fluctuation values. 

U2) For the first time step, the computation domain is solved as 
steady flow with the corresponding unsteady inlet velocity 
value. This is done because no previous results are available 
for the time dependent calculations and also it will provide 
initial guesses of the velocities, temperature, and pressure for 
the next time step. 

U3) The time step counter "t" advances by one and the 
computation domain is solved by the steps (S2) to (S10) from 
the steady flow procedure except the governing equations 
contain the time dependent, unsteady terms. The sinusoidal 
time variation of the inlet velocity profile is used in these 
calculations. 

U4) If the inlet flow reverses its direction on the next time step, 
the solutions are saved in a reverse order matter which 
simulates "flipping" the channel around in a physical sense. 


& /> 


Otherwise, the present solution is normally saved for the next 
time step calculations. 

U5) Steps (U3) and (U4) are repeated until the total number of 
time steps are completed. The results are printed for the last 
oscillating cycle. 




CHAPTER IV 


DISCUSSION OF RESULTS 


The presentation of the results are divided into three sections. The 
first section compares and validates the present steady flow solutions with 
existing literature on incompressible, thermally expandable and compressible 
fluid conditions. These variable fluid property validations determine and 
evaluate the accuracy of the present computer code prior to the 
examination of the pulsating flows and oscillating flows discussed in 
sections two and three, respectively. For reference purposes, the 
abbreviations and the numbers denoted in square brackets coincide with 
the flow runs performed in this chapter, and the detailed input parameters 
are summarized in Appendix-A. 


£ // 


Code Validation 


Steady Incompressible - Hydrodynamicallv Developing Flows : 

Based on Chiu’s [1984] thesis, the present computer code is first 
validated for its ability to produce accurate steady incompressible flow 
solutions prior to any code modifications on handling thermally expandable, 
compressible properties and unsteady conditions. In Figures 4.1.1a and 
4.1.1b, the residual mass source is plotted with respect to the number of 
"Global" iterations for the Reynolds numbers of 40 and 100, respectively. 
The figures reveal the convergence characteristic behaviors for steady 
incompressible flows at four axial locations from the area of flow 
development to the fully-developed downstream region. Having the over- 
relaxation factor for the pressure correction terms set at 1.5, the present 
code performed 53 "Global" iterations to reach the 0.2% mass flow 
convergence criterion at the two low Reynolds number flows considered. 
As for the high Reynolds number flows (Re=2000 or greater), substantial 
reduction in the required number of iterations (20 or less) was observed. 

The gradual velocity evolvements from the uniform inlet profile to 
the fully-developed, parabolic form are illustrated in Figures 4.1.2a and 
4.1.2b. For both of the Reynolds numbers (40 and 100), the fully- 
developed streamwise velocity reached a maximum normalized value of 1.5 
at the channel's center-line which meets the theoretical value for steady 
incompressible flows between parallel plates. Notice also that the near wall 
over-shoot velocities in the initial developing profiles correspond to 
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Figure 4.1.1: Convergence characteristics for incompressible, steady flow at: 

a) Re =40 [SI-1]. 

b) ReJ=100 [SI-2]. 
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through parallel plates at: 

a) Re =40 [SI-1]. 

b) Re c =100 [SI-2]. 
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Brandt’s [1966] findings and to those who have solved the Navier-Stokes 
equations. The over-shoot velocity comparison with data from Brandt are 
listed in Table 4.1.1; good agreements on the location and on the maxima 
value of the over-shoot velocity were established. 

Listed in Table 4.1.2 are center-line velocity values situated at 
various axial distance from the channel's entrance. Two Reynolds numbers 
(40 and 1000) were examined and compared with data from Brandt; the 
agreements of the results were well below half a percent difference. 

An exercise to determine the steady hydrodynamic entrance length 
was performed for three Reynolds numbers: 40, 100 and 1000. The 
hydrodynamic entrance length is defined as the required distance from the 
channel's inlet for U-velocity to reach 99 percent of the fully- developed 
value of 1.5. The outcomes are listed and compared with other 
literature in Table 4.1.3. 

The boundary layer theory results from Schlichting [1960] estimated 
shorter hydrodynamic entrance lengths than others who resolved the 
Navier-Stokes equations. At a Reynolds number of 40, Schlichting's 
prediction was 28.9% less than the present result; while at the higher 
Reynolds number of 1000, the difference decreased to 17.2% less than the 
present result. From this comparison, it is clear that the boundary layer 
theory provides satisfactory approximations of the Navier-Stokes equations 
only when the Re 0 is high. 

The present hydrodynamic entrance length results are in good 
agreement with data from Brandt [1966], Morihara [1973] and Narang [1976]. 



Table 4.1.1: Location and Amplitude of the maxima of the steady, developing 

velocity profiles compared with data from Brandt [1966] at different 
Re fl and axial locations [SI-1, SI-5]. 



Table 4.1.2: Center-line velocity comparison of steady, developing flow between 

parallel plates with data from Brandt [1966] at different Re. and 
axial locations [SI-1, SI-5]. 0 




U (at y/H=.5) 

Re„ 

x/H 

Brandt 

[1966] 

Present 

Diff. 


0.1 

1.022 

1.021 

- 0 . 1 % 

40 

0.2 

1.085 

1.079 

- 0 . 6 % 

0.3 

1.169 

1.166 

- 0 . 3 % 


0.4 

! 

1.254 

1.250 

- 0 . 3 % 


0.1 

1.003 

1.003 

0 . 0 % 

1000 

1.0 

1.100 

1.096 

- 0 . 3 % 

3.0 

1.211 

1.213 

0 . 2 % 


5.0 

1.288 

1.290 

0 . 2 % 


70 





















































































Table 4.1.3: Comparison of the steady hydrodynamic entrance length for flows 

in parallel plates with data from Schlichting [1960], Brandt [1966], 
Morihara [1973], Chen [1973] and Narang [1976]. [SI-1, SI-2, SI-5] 


Re, 

Hydrodynamic Entrance Length. x/H 
( at U=.99 *U -m ) 

Present 

Schlichting 

[1960] 

Brandt 

[1966] 

Morihara 

[1973] 

Chen 

[1973] 

Narang 

[1976] 

40 

1.125 

0.8 

1.13 

1.118 

1.50 

1.11 

100 

2.630 

2.0 

1 

1 

1 

1 

— 

2.91 

2.72 

1000 

24.15 

20.0 

22.4 

— 

l 

l 

1 

i 

24.42 


Methods used by the earlier investigators: 

Schlichting : Boundary Layer Theory 
Brandt : Full Navier-Stokes Equation 
Morihara : Numerical Solution, Quasilinearized 
Chen : Momentum Integral Method 

Narang : Linearizing the Inertia Terms 


Table 4.1.4: Steady mean nusselt number for equal and constant walls' 

temperature at Re =100 and Pr=.72 compared with data from Hwang 
[1973]. [SI-2] 


| 

Re. 

x/H 

Mean Nil* 

Hwang [1973] 

Present 

Diff. 

100 

.500 

14.90 

16.74 

12.3% 

1.25 

10.96 

11.72 

6.9% 

2.13 

9.593 

10.03 

4.6% 
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As for Chen [1973], who solved the Navier-Stokes equations by the 
momentum integral method, he predicted longer hydrodynamic entrance 
lengths than the present work (at Re 0 =40 which was 33% longer, and at 
Re 0 =lOO which was 10.6% longer). 


Steady Incompressible - Simultaneously Developing Flows : 

In this segment the heat transfer analysis in the channel, provided 
by the constant wall temperature is examined with both velocity and 
temperature profiles developing simultaneously. 

The first validation case has a constant temperature for the upper 
and lower walls while the flow is at Re Q =100 and the fluid has a Pr=.72. 
The temperature evolvement from the inlet uniform profile to the fully- 
developed parabolic shape is illustrated in Figure 4.1.3a, while Figure 
4.1.3b plots the present heat transfer coefficient results along with the 
available numerical data from Hwang [1973]. The Nusselt numbers are in 
good agreement between the two predictions with a maximum difference of 
12% after the channel's length x/H=.5 (see Table 4.1.4). Moreover, the 
present local Nusselt number curve approaches the analytical fully- 
developed value of 7.54 downstream. 

The second heat transfer validation run has constant and different 
wall temperatures, where the upper plate heats the fluid while the lower 
plate cools it. Figure 4.1.4a shows the results of the developing 
temperature profiles along 5 axial locations from the entrance to the fully- 
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developed downstream region. The developing temperature profiles began 
uniformly and then as the input heat flux from the upper wall balanced 
with the heat flux ejected into the lower wall, the temperature profile 
transformed into a fully-developed diagonal linear line with the end-points 
equal to the upper and lower wall temperatures. In Figure 4.1.4b, the local 
Nusselt number curve is plotted and the curve approaches the fully- 
developed theoretical value of 4.0. 


Steady Thermally Expandable Flows : 

For the thermally expandable flows, the work done by Schade and 
McEligot [1971] is used as the comparison reference. In their paper, the 
predictions of the fluid properties were determined by the ideal gas law, 
the viscosity power law and the conductivity power law. These equations 
are defined in Chapter 2 (eqs. 2.9, 2.10 and 2.12) with the a and 0 values 
listed in Table 2.1 for air. 

Two Reynolds numbers, 144 and 2000, were chosen for this extensive 
validation study, and the fluid was assumed to have the properties of air 
and four T W /T Q ratios (.5 cooling, 1.0 incompressible, 2.0 and 5.0 heating) 
were examined for each Re 0 . The resulting comparisons are shown in 
Figures 4.1.5 to 4.1.7. 

In Figures 4.1.5a (Re Q =144) and 4.1.5b (Re 0 =2OOO), the local apparent 
friction is compared with data from Schade and McEligot at T w /T 0 =1.0 
(incompressible), 0.5 and 2.0 (thermally expandable). The present 







Figure A. 1.5: Local apparent friction comparison of steady flow in the channel 

with data from Schade & McEligot [1971] at different wall/inlet 
temperature ratios, T„/T # »1.0 (incompressible), 0.5 and 2.0 
(thermally expandable): 

a) Re =144 [SI-4, ST-1, ST-2]. 

b) Re =2000 [SI-5, ST-4, ST-9]. 
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Ratio of wall and bulk-mean temperature comparison in the 
channel with data from Schade & McEligot [1971] at different 
wall/inlet temperature ratios, T^/T^O.5, 2.0, 5.0 (thermally 
expandable): 
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.1.7: Local Nusselt number comparison of steady flow in the channel 

with data from Schade & McEligot [1971] at different wall/inlet 
temperature ratios, Ty/T 9 =1.0 (incompressible), 0.5 and 2.0 
(thermally expandable): 

a) Re =144 [SI-4, ST-1, ST-2]. 

b) Re°=2000 [SI-5, ST-4, ST-9]. 




predictions suggest an increase in friction as the fluid is heated while a 
decrease is observed when cooled; this observation corresponds to Schade 
and McEligot’ s findings. This friction behavior is directly related to the 
increase in viscosity as the temperature of the gas rises and vice versa. 
Such variation in friction can only be recorded when the fluid properties 
are allowed to vary; whereas for incompressible flows, a single friction 
curve would have been generated regardless of the wall/inlet temperature 
ratio. A quantitative comparison revealed that the present friction 
predictions compared well with data from Schade and McEligot [1971] 
especially at the higher Re 0 (2000), (see Table 4.1.5). 

The wall and bulk-mean temperature ratios versus the normalized 
axial locations are plotted in Figures 4.1.6a (Re 0 =144) and 4.1.6b (Re 0 =2000). 
Again, the comparisons with data from Schade and McEligot are in excellent 
agreement, even at high heating of T y /T 0 =5.O. 

The local Nusselt number comparison at wall/inlet temperature ratios 
of 1.0 (incompressible), 0.5 and 2.0 (thermally expandable) are shown in 
Figures 4.1.7a (Re 0 =144) and 4.1.7b (Re 0 =2OOO). Notice that the Nusselt 
number curve increases as a restilt of a higher wall/inlet temperature ratio 
and the opposite effect is shown for cooling. This behavior is caused by 
the fluid's conductivity increase with the temperature rise which offers 
less thermal resistance against heat flow into the fluid. Quantitative 
comparison with Schade and McEligot's data resulted in an agreement 
within a 5.0% difference (see Table 4.1.5). 


Table 4.1.5: Comparison of the predictions for steady, thermally expandable flow 

between the present and the data from Schade & McEligot [1971]. 
[ST-1, ST-2, ST-4, ST-9] 


Re. 

T./T. 

B 

* Re, / 24.0 

T. / T, 

Nu , | 

Schade 

•• 

Present 

Diff. 

(%) 

Schade 

•• 

Present 

Diff. 

(%) 

Schade 

• • 

Present 

Diff. 

(%) 

144 

1 

.01 

1.48 

1.68 

13.5 

.541 

.548 

1.3 

10.0 

10.3 

3.0 

.05 

.872 

.911 

13 

.609 

.621 

El 

7.12 

7.42 

■a 

.10 

.806 

.874 

8.4 

.676 

.690 

KB 

7.05 

7.40 

d 

.20 

.841 

.944 

12.2 

.776 

.798 

2.8 

7.21 

7.49 

3.9 

1 

.01 

2.49 

2.69 

8.0 

1.70 

1.65 

- 2.9 

11.2 

11.0 

- 1.8 

.05 

1.40 

1.31 

- 6.4 

1.37 

1.34 

- 2.2 

8.11 

8.01 

- 1.2 

.10 

1.17 

1.12 

EB 

1.19 

1.17 

- 1.7 

7.81 

7.80 

- 0.1 

.20 

1.05 

1.01 

- 3.8 

1.06 

* 1.05 

- 0.9 

7.62 

7.68 

0.8 

2000 

B 

.01 

1.48 

1.54 

13 

.541 

.545 


10.0 

10.3 

^^6 

.02 

1.14 

1.18 

3.5 

.561 

.568 

1.2 

8.22 

8.55 

■a 

.05 

.872 

.883 

1.3 

.609 

.619 

1.6 

7.12 

7.43 

■a 


.01 

2.49 

2.45 

- 1.6 

1.70 

1.68 

- 1.2 

11.2 

10.8 

- 3.6 

.02 

1.90 

1.81 

ESI 

1.58 

1.56 

— 1.3 

9.27 

9.09 

- 1.9 

.05 

1.40 

1.35 

- 3.6 

1.37 

1.36 


8.11 

8.06 

- 0.6 


Note: ** - Schade It McEligot [l97l] used the Power Laws for properties variations. 








































































































































Steady Compressible Flows : 


When dealing with the steady high velocity flow, previous 
investigators mainly studied external supersonic flow or subsonic flow 
over a flat plate or a cross section of an air foil. To the author's 
knowledge, no internal subsonic compressible flow through parallel plates 
was published. Hence, the results in this segment are new findings for 
internal compressible flows. 

Ensuring the code generates accurate compressible solutions, 
validation runs were established at the fixed Re Q =2000 and M=0.25 tested 
under four different wall/inlet temperature ratios. Note that when dealing 
with compressible flow, the stagnation enthalpy or the stagnation 
temperature (since C p is assumed constant) is taken as the thermal energy 
reference of the fluid. By setting the inlet fluid stagnation temperature 
equal to 300K, the four different wall temperatures (note: the static wall 
temperature is equivalent to the stagnation wall temperature due to zero- 
velocity at the wall) were set constant at 294K and 298K (fluid cooling, 
cases A and B, respectively), 300K (no heat transfer, case C) and 302K 
(fluid heating, case D). The static and stagnation temperature profiles are 
presented in Figures 4.1.8a to 4.1.8d where each figure represents 
solutions at a specific streamwise location: a) x/H=1.10, b) x/H=5.16, c) 
x/H=19.9 and d) x/H=59.7. 

The present computer code handled the static temperature as one of 
the primitive variables (see eq. 2.5); therefore, the direct results from the 
energy equation are static temperature solutions which are plotted as solid 
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lines in figure 4.1.8a to 4.1.8d. Notice for all four T v /t q ratios examined, 
the static temperature results (solid lines) gradually decreased in value 
with respect to the axial length. By intuitive observation, one would 
conclude that there was fluid energy loss into the wall at all four of the 
wall/inlet temperature ratios studied. However, by close examination of the 
temperature gradients at the wall suggested otherwise; gradients for cases 
A and B indicate negative heat flux to the fluid (cooling), case D revealed 
a positive heat flux to the fluid (heating), while case C had a zero wall 
temperature gradient hinting no heat was transferred. The same 
conclusions can also be derived by viewing the stagnation temperature 
curves represented as the dashed lines in the plots. 

As explained in Chapter two, the difference between the static and 
stagnation temperature curves is that the static temperature measures only 
the thermal energy in the fluid while the stagnation temperature accounts 
for the total energy (thermal and kinetic energy). Whereas, in the 
incompressible and thermally expandable fluid assumptions, the thermal 
energy outweighs the kinetic energy. The two types of energy in the high 
velocity compressible flow are comparable. For this very reason, the 
decrease of static temperature with respect to the axial length was caused 
by the conversion of thermal energy into kinetic energy, and if the static 
temperature was examined alone, it could have led to a false conclusion on 
the heat flux direction. This validation analysis agrees with the similar 
qualitative findings for high speed external flow over a flat plate [Kays 
and Crawford, 1980]. 
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The main objective of this thesis is to explore the fluid 
compressibility effects on internal flows as compared with the 
incompressible and thermally expandable assumptions. The following study 
examines the behavior of the apparent friction factor and the heat transfer 
coefficient under various Mach number flows. The fluid air was selected 
for the investigation with the Re Q fixed at 2000. Two sets of steady flow 
with heat transfer were performed, one set on fluid heating (T^/r^l.5) and 
one set on fluid cooling (T^/r^.67). Five fluid property models were 
considered: a) incompressible; b) thermally expandable; c) compressible, 
M=.05; d) M=.10 and e) M=.25. 

The predictions of the apparent friction are plotted in Figure 4.1.9a 
(for the heating case) and Figure 4.1.9b (for the cooling case). As 

expected, substantial differences exist between the incompressible and 
thermally expandable predictions. For the T^r fl =1.5 heating case, the 
thermally expandable run predicted around 12.5% higher in friction than 
the incompressible results throughout the channel; while for the T tf /r 0 =.67 
cooling case, about 15.5% lower was found (see Table 4.1.6). As noted in 
the previous thermally expandable flow discussion, this behavior is 
primarily caused by the temperature dependence of the fluid viscosity (eq. 
2 . 11 ). 

Since thermally expandable flow is actually compressible flow with the 
Mach number assumed as zero, the friction factor comparison with 
compressible flows at low Mach numbers shows marginal differences. For 
both fluid heating and cooling studies, thermally expandable results differ 
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Figure 4.1.10: 


Comparison of steady local Nusselt number in the channel at 
different Mach numbers, M=0.0 (incompressible), 0.0 (thermally 
expandable), .05, .10 and .25 (compressible): 

a) Vr =1.5 [SI-8, ST-8, SC-8, SC-9, SC-10]. 

b) TVtI=.67 [SI-6, ST-5, SC-5, SC-6, SC-7]. 
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Table 4.1.6: Steady apparent friction comparison between incompressible, 

thermally expandable, and compressible air flow for Re =2000 and 
different Mach numbers [SI-6, SI-8, ST-5, ST-8, SC-5° to SC-10]. 


T./t. 

x.* 

F.„*Re 1 /24.0 

Diff. from T.E. (%) 

Inc. 


Mach No. 

Inc. 

Mach No. 

.05 

.10 

.25 

.05 

.10 

.25 

1.5 

.02 

1.449 

1.662 

1.665 

1.675 

1.762 

-12.8 

0.2 

0.8 

6.0 

.04 

1.149 

. 

1.323 

1.327 

1.339 

1.449 

-13.2 

0.3 

1.2 

9.2 

.06 

1.051 

1.199 

1.203 

1.215 

1.344 

-12.3 

0.3 

1.3 

12.1 

.08 

1.019 

1.145 

1.149 

1.161 

1.309 

-11.0 

0.3 

1.4 

14.3 

.667 

.02 

1.449 

1.258 

1.259 

1.262 

1.288 

15.2 

0.1 

0.3 

2.4 

.04 

DS3 

0.993 

0.994 

0.998 

1.027 

15.7 

0.1 

0.5 

EZl! 

.06 

1:051 

0.906 

0.908 

0.911 

0.942 

16.0 

0.2 

0.6 

4.0 

.00 

1.019 

0.881 

0.882 

0 886 

0.919 

15.7 

0.1 

0 6 

4.3 


Table 4.1.7: Steady Nusselt number comparison between incompressible, 

thermally expandable, and compressible air flow for Re =2000 and 
different Mach numbers [SI-6, SI-8, ST-5, ST-8, SC-5 to SC-10]. 


T./t. 

x 0 * 

Nu b 

Diff. from T.E. (%) 

Inc. 

m 

Mach No. 

Inc. 

Mach No. 

.05 

.10 

.25 

.05 

.10 

.25 

1.5 

.02 

0.819 

9.000 

9.017 

9.044 

9.273 

-2.1 

0.1 

0.4 

2.9 

.04 

7.960 

8.183 

8.194 

8.228 

0.533 

-2.6 

0.1 

0.5 

4.3 

.06 

7.820 

8.027 

8.041 

8.083 

0.505 

-2.6 

0.2 


6 0 

.08 

7.798 

7.977 

7.994 

8.046 

0.620 

-2.2 

0.2 

0.9 

8.1 

.667 

.02 

8.819 

8.563 

8.560 

8.552 

0.495 

3.0 

-0.1 


-0.8 

.04 

7.968 

7.680 

7.676 

7.666 

7.593 

3.8 

-0.1 


-1.1 

.06 

7.820 

7.520 

7.516 

7.504 

7.411 

4.0 

-0.1 


-1.4 

.08 

7.798 

7.510 

7.505 

7.491 

7.382 

3.8 

-0.1 


-1.7 


Note: 

Inc. - Incompressible 

T.E. - Thermally Expandable 

Mach No. - Compressible flow at the given Mach number 













































































































































































less than 1.5% when compared with compressible flows at Mach numbers 
below .10 . However as the Mach number increases to .25, the compressible 
friction factors surpass the thermally expandable predictions regardless of 
the wall/inlet temperature ratios; Table 4.1.6 discloses the differences that 
were as high as 14.3% (at X 0 + =.O8 and T^r^l.5). 

Next, the Nusselt number predictions are plotted in Figure 4.1.10a 
(heating) and Figure 4.1.10b (cooling). Similar behaviors as found in the 
friction factor apply to the Nusselt number except the differences between 
compressible and thermally expandable results are lower in percentage 
magnitude (see Table 4.1.7). 

Via this steady compressible flow study, the following observation 
was noted. For low Mach number flows (below M=.10), the thermally 
expandable fluid assumption predicts excellent results as compared to the 
actual compressible flow. However, as the Mach number of the flow 
increases (above M=.25), the thermally expandable predictions will no longer 
be valid to model such high-velocity compressible flow. 


Pulsating Flows 

With the preliminary steady flow investigation conducted and the 
code validated, the next research stage implements the time-dependent 
terms into the governing equations. In this section, the periodic pulsating 
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unsteady flow is examined in detail. In part, the present pulsating 
incompressible and thermally expandable results were presented and 
published at the fourth International Symposium on Transport Phenomena 
in Heat and Mass Transfer at Sydney, Australia [Ibrahim and Kwan, 1990]. 


Pulsating Incompressible - Hvdrodynamically Developing Flows : 

The code modifications for unsteady flow are first checked with the 
incompressible pulsating results from Siegel and Perlmutter [1962]. The 
pulsating flow boundary conditions described in Chapter two are similar to 
the boundary conditions employed by Siegel and Perlmutter, except they 
applied the sinusoidally pulsating pressure gradients at the channel's inlet 
rather than the pulsating velocity (eq. 2.14). Ensuring the two inlet 
velocity matches, Siegel and Perlmutter's equations were used to derive the 
present imposed velocity fluctuation amplitude, ip. The ip's are listed in 
Table 4.2.1 for various Valensi numbers with the imposed pressure gradient 
fluctuation amplitude equal to one. 

The predicted phase angle between the axial pressure gradient and 
the mean velocity are compared between the present work versus the 
closed form incompressible solution given by Siegel and Perlmutter [1962] 
in Table 4.2.2. For the low and high Valensi numbers considered, the 
comparison shows excellent agreement. 

At the fixed Re re j=2000, the Figures 4.2.1a (o*=0.08) and 4.2.1b 
(o>*=32.0) illustrate the fully developed, normalized fluctuating velocity 
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Table 4.2.1: Parameter f used in determining the inlet pulsating velocity 

conditions (eq. 2.14) for different to match the conditions 
used by Siegel & Perlmutter [1962] [PI-1 to PI-4]. 


cj' 


0.08 

1.0000 

0.0 

0.7775 

32.0 

0.2984 

200.0 

0.0543 


Table 4.2.2: Comparison between present work and data from Siegel & 

Perlmutter [1962] on the phase angle difference between cross- 
sectional average velocity and axial pressure gradient 
(incompressible pulsating flow, fully developed) [PI-1 to PI-4]. 


w* 

Phase Angle Difference 
(degrees) 

Siegel & Perlmutter 
[1962] 

Present 

0.08 

0.0 

0.0 

8.00 

38.5 

38.7 

32.0 

70.5 

71.0 

200.0 

83.5 

83.1 


9t 


















y/H 

Figure 4.2.1: Cross sectional fluctuating velocity profiles compared between 

present work and data from Siegel & Perlmutter [1962] (fully 
developed, incompressible pulsating flow, Re |eja =2000): 

a) tt*=0.08, *=1.00 [PI-1]. 

b) o>*=32.0, *=0.298 [PI-3]. 




=200.0 f =0.0543 







versus the normalized distance from the wall. The figures are comparisons 
between the present work and data from Siegel and Perlmutter; excellent 
agreement was obtained. 

A study of the developing velocity profile for pulsating flow was 
carried out for Re re j=2000 and <i>*=200. Figures 4.2.2a to 4.2.2c show the 
predicted normalized fluctuating velocity versus normalized distance from 
the wall at three different axial locations. An overshoot has been noticed 
in the velocity profile near the duct inlet which dissipates downstream. 
This observation is similar qualitatively to Creff’s et.al. [1983,1985] findings 
for pulsating flow in circular tubes. 


Pulsating Incompressible - Thermally Developing Flows : 

(Slug Flow Approximation) 

Figures 4.2.3a and 4.2.3b show the predicted Nusselt number versus 
the normalized axial distance comparing the present work with data from 
Siegel and Perlmutter [1962] (Slug flow approximation). The Valensi 
numbers considered are 0.08 (Figure 4.2.3a) and 32.0 (Figure 4.2.3b). It 
should be noted that Figure 4.2.3a shows the instantaneous Nusselt number, 
while Figure 4.2.3b shows the fluctuating Nusselt number (instantaneous 
minus steady state value). The figures clearly illustrate that the 
agreement between the present work and the closed form solution for 
thermally developing with slug flow approximation is excellent. It can be 
seen from the figures that, at a low Valensi number, the Nusselt number 
decreases monotonically as x increases (at all values of cat - see Figure 
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Figure 4.2.3: Nusselt number variation in the channel compared between 

present work and data from Siegel & Perlmutter [1962] 
(incompressible with pulsating slug flow approximation, 
Re ieair 2000 ’ Pr=0-72): 

a) Instantaneous Nusselt number, o> =0.08, ^=1.00 [PI-1]. 

b) Fluctuating component of the Nusselt number, <■> =32.0, 

*=0.298 [PI-3]. 




4.2.3a); accordingly, the quasi-steady approximation can be applied. As the 
Valensi number increases the fluctuating Nusselt number takes a complex 
shape in its variation with x or time. The quasi-steady approximation can 
only be applied close to the duct inlet. 


Pulsating Simultaneously Developing Flows (Variable Fluid Properties) : 

In this segment, new results will be presented for thermally and 
hydrodynamically developing flow for incompressible, thermally expandable 
and compressible fluid assumptions. The effects of variable fluid 
properties and fluctuating frequencies (Valensi number) on the pulsating 
flow results are discussed. The following parameters are fixed for the rest 
of the investigation in this pulsating flow solution section: a) Re fef =2000, b) 
air is the fluid, and c) i|r=1.0 (50% of the mean velocity value). 

Figures 4.2.4a to 4.2.4c show the mass flux variation in the channel 
at different cat. Figure 4.2.4a illustrates the results for incompressible flow 
at any g>* value. Notice the inlet mass flux is always equivalent to the exit 
mass flux at any instant of o>t. This behavior is due to the constant 
density restriction applied by the incompressible fluid assumption; the same 
conclusion was observed for compressible flow at a low Valensi number 
(<■>*=0.08). But as the Valensi number increases, the instantaneous mass flux 
at the two ends of the duct is no longer equal for variable property flows. 
The viscous damping of the fluid reduces the amplitude of the unsteady 
mass flux fluctuation downstream (see Figures 4.2.4b and 4.2.4c). 
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Mass flux variation in the channel with different fluid 
property assumptions (incompressible, thermally expandable 
and compressible), (Re„ =2000, T /t =1.2, Pr=.72, *=1.00): 











Accordingly, Table 4.2.3 presents the Amplitude of the exit mass flux 
fluctuation for the three a* (.08, 32.0 and 100.0) examined with the inlet 
fluctuation amplitude set equal to 1.0 . The following pulsating flow 
observations were gathered from Figures 4.2.4a to 4.2.4c and Table 4.2.3: 

a) For incompressible flows, the instantaneous mass flux is 
identical at any axial location which means any changes or 
disturbance in the inlet mass flux is immediately transferred 
throughout the channel. This is true at any Valensi number for 
pulsating incompressible flow. 

b) For thermally expandable and compressible flows, the 
(. fluctuating mass flux signal at the inlet is delayed downstream by 

the fluid density changes. Also the amplitude of mass flux 
fluctuation decreases with respect to the axial distance caused by 
the fluid viscous damping. These results suggest that if the channel 
was sufficiently long, the damping effect of the fluid will eventually 
eliminate the fluctuating component of the pulsating flow; hence, the 
flow will become steady-state downstream. 

c) Results in Table 4.2.3 suggest the higher the Valensi number 
(dimensionless frequency) of the pulsating flow, the faster the 
fluctuating component diminishes. This phenomenon is similar to the 
behavior of sound wave, where lower frequency waves can travel 
further in distance than higher frequency waves at the same 
amplitude of fluctuation. 
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Table 4.2.3: Amplitude of the fluctuating mass flux at the channel's exit for 

different <■>* (thermally expandable pulsating flow, L/H=70.0, 
Re lean =20° 0 , Pr=.72, T/T # =1.2, *=1.00) [PT-1, PT-3, PT-5]. 


CJ- 

( G eirtt / G„ f ) m „ 


Amplitude of the 
Exit Mass Flux 
Fluctuation 

.08 

1.500 

0.500 

1.000 

32.0 

1,489 

0.527 

0.962 

100.0 

1.439 

0.555 

0.884 


I 


Table 4.2.4: Comparison between the time-averaged pulsating friction with 

the steady friction at different wall/inlet temperature ratios 
{Re J 83 J =2000 » Pra - 72 > C pI -5 to PI-7, PT-1 to PT-6]. 


T,/T. 

! 

Steady 

Time-Averaged Local F* 

Diff. from the Steady F* 
(%) 



Local 








■1 


u*=.08 

u*=32 

u*=100 

e 

• 

n 

o 

CD 

cj*=32 

cj*=100 


.005 

2.619 

2 589 

2.590 

2.599 

wbsm 

-l.l 

-0.8 

1.0 

(lac.) 

.010 

1.914 

1.890 

1.891 

1.896 

flfl 

-1.2 

-0.8 

.020 

1.449 

1.434 

1.433 

1.424 

■ 

-1.1 

-1.7 

.040 

1.149 

1.152 

1.142 

' 1.098 


-0.6 

-4.4 


.080 

1.019 

1.031 

1.012 

0.936 

1.2 

-0.7 

-8.2 


.005 

2.785 

2.753 

2.753 

2.754 

-1.1 

-1.1 

-1.1 

1.2 

(T.E.) 

.010 

2.037 

2.012 

2.011 

1.996 

-1.2 

-1.3 

-2.0 


1.543 

1.521 

1.515 

1.470 

— 1.4 

-1.0 

-4.7 

'tT 9 

1.225 

1.208 

1.191 

1.110 

-1.4 

-2.8 

-9.4 


ESI 

1.078 

1.064 

1.037 

0.962 

-1.3 

-3.8 

-10.8 


.005 

2.931 

2.894 

2.693 

2.886 

-1.3 

-1.3 

-1.5 

1.4 

(T.E.) 

.010 

2.143 

2.116 

2.112 

2.083 

-1.3 

-1.4 

-2.8 

2 

1.624 

1.596 

1.587 

1.516 

-1.8 

-2.3 

-6.5 

isp! j 

1.291 

1.260 

1.238 

1.137 

-2.4 

-4.1 

-11.9 



1.123 

1.095 

1.065 

0.992 

-2.5 

-5.2 

-11.7 


Note: F* - * Re, / 24 

Inc. - Incompressible 
T.E. - Thermally Expandable 


99 













































Figures 4.2.5a to 4.2.5c present the predicted local time-averaged 
apparent friction versus the normalized axial distance with ti>*=0.08 (Figure 
4.2.5a), 32.0 (Figure 4.2.5b) and 100.0 (Figure 4.2.5c). The pulsating results 
are plotted with the steady results as reference. In addition, Table 4.2.4 
shows the numerical comparison between the steady state and the local 
time-averaged apparent friction under the incompressible and thermally 
expandable models. For a>*=0.08 and 32.0, the difference between the steady 
state friction and the local time-averaged friction is within 5% for the 
three T^/tp (1.0, 1.2 and 1.4) ratios. The results lead to the conclusion 
that these low frequency flows are still close to quasi-steady. However, 
at a Valensi number equal to 100.0, the difference between the steady state 
friction and the local time-averaged friction increases to about 12% 
downstream. Only the near entrance region can be modeled by the quasi- 
steady assumption. 

Comparison of the local time-averaged friction between thermally 
expandable and compressible models is presented in Figures 4.2.6a to 4.2.6c. 
At the lower frequencies (<o*=0.08 and 32.0), the difference is minimal (see 
Table 4.2.5), but at the higher frequency ( o>*= 100.0) the difference between 
the two fluid models reached as high as 7.6% downstream, 

Figures 4.2.7a to 4.2.7c show the predicted normalized mean fluid 
temperature versus time at different axial locations and <i>*=0.08 (Figure 
4.2.7a), 32.0 (Figure 4.2.7b), 100.0 (Figure 4.2.7c). As shown in the figures, 
the temperature fluctuation increases as the fluid moves downstream, while 
the amplitude of these temperature fluctuations decreases as the Valensi 
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Figure 4.2.5 (cont.): 

c) o*= 100.0 [PI-7, PT-5, PT-6]. 










Table 4.2.5: Comparison of the time-averaged pulsating friction between 

predictions from thermally expandable and compressible flows 
(Re^^OOO, air) [PT-1, PT-3, PT-5, PC-1 to PC-3]. 


0>* 

T./t. 

Mach 

no. 

Y * 

A r*( 

Time-Averaged Local F* 

Diff. between Comp, 
and T.E. (%) 

T.E. 

Comp. 

.08 

1.2 

.05 

.005 

2.753 

2.756 

0.1 

.010 

2.012 

2.015 

0.1 

.020 

1.523 

1.526 

0.2 

.040 

1.210 

1.212 

0.2 

.080 

1.066 

1.068 

0.2 

32 

1.2 

.05 

.005 

2.753 

2.750 

0.2 

.010 

2.011 

1.999 

-0.6 

.020 

1.515 

1.496 

-1.3 

.040 

1.191 

1.178 

-1.1 

.080 

1.037 

1.015 

-2.1 

100 

1.2 

.02 

.005 

2.754 

2.755 

0.0 

.010 

1.996 

1.996 

0.0 

.020 

1.470 

1 462 

-0.5 

.040 

1.110 

1.077 

-3.0 

.080 

0.962 

0.889 

-7.6 


Note: F* = • Re, / 24 

T.E. - Thermally Expandable 
Comp. - Compressible 
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number increases. Moreover, for the same first order harmonic flow field, 
the temperature field contains only low order harmonics at low Valensi 
numbers (see Figure 4.2.7a); while at higher Valensi numbers (see Figures 
4.2.7b and 4.2.7c), distortions of the temperature fluctuations suggest 
higher order harmonics are present. The trend of the temperature field 
is similar for both incompressible and thermally expandable models, except 
the incompressible model under-states the temperature. Table 4.2.6 lists 
the numerical comparison of Figures 4.2.7a to 4.2.7c. 

Figure 4.2.8 shows similar plots of mean temperature fluctuation as 
presented in Figure 4.2.7 except the predictions are from thermally 
expandable and compressible fluid models. The trend of the temperature 
field is identical for both models with minimal difference. 

The effect of the Valensi number towards the heat transfer of the 
pulsating flow is examined and presented in Figures 4.2.9a (o)*=0.08), 4.2.9b . 
(<i>*=0.08) and 4.2.9c (<i>*=0.08). The figures plot the local time-averaged 
Nusselt number together with the corresponding steady Nusselt number 
results at different wall/inlet temperature ratios. The comparison of the 
time-averaged results versus the steady Nusselt number results shows 
increased differences as the Valensi number increases, however, even at 
the high *=100.0, the maximum difference downstream is only 1.2% (see 
Table 4.2.7). Thus, the steady heat transfer predictions are good 
estimations of the time-averaged pulsating heat transfer results even at 
high Valensi numbers. 

Illustrated in Figures 4.2.10a to 4.2.10c are heat transfer solutions 
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Table 4.2.6: Periodic location and amplitude of the maxima cross-sectional 

mean fluid temperature [PI-5 to PI-7, PT-2, PT-4, PT-6]. 


«• 

T./T. 

x/H 

Maximum 

Velocity Phase Angle, ut 
at the Maximum & m ^ 
(degrees) 

Incomp. 

T.E. 

Incomp. 

T.E. 

.08 

1.4 

1.1 

1.060 

1.060 

270 

270 

5.2 

1.110 

1.130 

270 

270 

20.0 

1.206 

1.221 

270 

270 

45.0 

1.287 

1.305 

270 

270 

60.0 

1.310 

1.336 

270 

270 

32 

1.4 

1.1 

1.060 

1.068 

270 

270 

5.2 

1.118 

1.130 

280 

275 

20.0 

1.203 

wmm 

300 

295 

45.0 

1.271 

1.293 

330 

320 

60.0 

1.296 

1.319 

340 

335 

100 

1.4 

1.1 

1.060 

1.060 

275 

275 

5.2 

1.117 

1.130 

295 

295 

20.0 

1.188 

1.200 

350 

345 

45.0 

1.236 

1.259 

45 

40 

60 0 

1254 

1.278 

75 

70 


Notes: Incomp, - Incompressible 

T.E. - Thermally Expandable 













































































.01 .03 .05 .07 .09 



Figure 4.2.9: Time-averaged Nusselt number in the channel at different 

wall/inlet temperature ratios: T^/T =1.0 (incompressible), 1.2 and 
1.4 (thermally expandable); (Re aB3n =2000, Pr=0.72, ^r=1.00): 

a) o»*=0.08 [PI-5, PT-1, PT-2]. 

b) <i>'=32.0 [PI-6, PT-3, PT-4], 








Table 4.2.7: Comparison between the time-averaged pulsating Nusselt number 

with the steady Nusselt number at different wall/inlet 
temperature ratios (Re IgaD =2000, Pr=.72) [PI-5 to PI-7, PT-1 to 

PT-6]. 


T./T. 

Y ♦ 

A r«f 

Steady 

Time-Averaged Local Nu k 

Diff. from the Steady Nu* 

Local Nu k 

<y*=. 08 

u*=32 

cj*=100 

«'=.08 

<u*=32 

cj*=100 


.02 

8 82 

8.83 

8.02 

8.84 

.11% 

0 . 0 % 

. 23 % 

1.0 

.04 

7.97 

8.02 

8.02 

7.99 

. 03 % 

. 63 % 

. 25 % 

(Inc.) 

.06 

7 02 

7.86 

7.06 

7.79 

. 51 % 

. 51 % 

-. 38 % 


.08 

7.80 

7.01 

7.02 

7.72 

. 13 % 

. 26 % 

- 1 . 2 % 

mm 

.02 

8.92 

8.91 

8.91 

8.92 

-.11% 

-.11% 

0 . 0 % 

IK 

.04 

0.08 

8.11 

8.11 

8.09 

. 37 % 

. 37 % 

. 12 % 

(T.E.) 

.08 

7.93 

7.94 

7.95 

7.91 

. 13 % 

. 25 % 

-. 25 % 


.08 

7 . B 9 

7.89 

7.90 

7.84 

0 . 0 % 

. 13 % 

-. 63 % 


.02 

0.98 

0.97 

8.97 

8.96 

-.11% 

-.11% 

1 

1.4 

.04 

8.15 

0.17 

0.17 

8.17 

. 25 % 

■ 


(T.E.) 

.06 

8.00 

7.99 

8.00 

7.99 

-.12% 


-. 12 % 


08 

7.95 

7.94 

7.95 

7.92 

-. 13 % 


-. 38 % 


Note: Inc. 
T.E. 


Incompressible 
Thermally Expandable 









































Table 4.2.8: Comparison of the time-averaged pulsating Nusselt number 

between predictions from thermally expandable and compressible 
flows (Re^ZOOO, air) [PT-1, PT-3, PT-5, PC-1 to PC-3]. 


u* 

T./r. 

Mach 

no. 

H 

Time-Averaged Local Nu, 

Diff. between Comp, 
and T.E. (%) 

T.E. 

Comp. 

.08 

1.2 

.05 

.010 

10.59 

10.59 

0.00 

.020 

8.912 

8.907 

0.06 

.040 

8.107 

8.099 

0.10 

.080 

7.889 

7.878 

0.14 

32 

1.2 

.05 

.010 

10.59 

10.55 

0.04 

.020 

8.911 

8.886 

0.28 

.040 

8.109 

8.091 

0.22 

.080 

7.901 

7.871 

0.38 

100 

1.2 

.02 

.010 

10.60 

10.58 

0.19 

.020 

8.922 

8.893 

0.33 

.040 

0.094 

8.042 

0.65 

.080 

7.839 

7.786 

0.68 


Notes: T.E. - Thermally Expandable 

Comp. - Compressible 


























































of the pulsating flow influenced by the effect of fluid compressibility 
under different pulsating frequencies. As listed in Table 4.2.8, similar 
Nusselt number predictions were found by the thermally expandable model 
and the compressible model (Mach number of .02 and .05). 

A couple of major observations are identified from this pulsating flow 
study. The fluid compressibility has a direct effect on the friction factor 
and there is a secondary effect on the heat transfer coefficient (see Table 
4.2.5 and 4.2.8). At a range of the Valensi number, the steady flow model 
provides good estimations for the time-averaged pulsating Nusselt number; 
whereas, for the time-averaged pulsating friction factor, close estimations 
by the steady flow model is limited to pulsating flow at the low Valensi 
number. 


Oscillating Flows 

The examination of the oscillating flow is divided into two segments. 
The first segment validates the code by comparing results with data from 
the existing literature on the incompressible model, then it leads into the 
investigation of compressibility effects on the friction and heat transfer 
coefficients. The second segment predicts friction and heat transfer 
results under the NASA Space Power Research Engine (SPRE) boundary 
conditions by the incompressible and compressible models. The comparison 
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of the two models will be presented and discussed. 


Oscillating Flow ( os =32.0 and Re m 3 a =2000) : 

Kurzweg [1985] developed an analytical equation to calculate the 
velocity profile existing in the central channel for incompressible, laminar, 
oscillating flow conditions. A validation incompressible run was made at 
(o=32.0 and R 6^=2000. The agreements between the present velocity profile 
and predictions from Kurzweg are excellent, as shown in Figure 4.3.1. The 
velocity profile pattern is symmetrical around the zero velocity which 
signifies that the flow is fully-developed. 

Figure 4.3.2 illustrates the pressure in the channel at various <ot. 
Besides the entry effect at both ends of the channel, the pressure 

fluctuation is symmetrical around the normalized pressure value of zero; 

J99o 

the same qualitative finding was observed by Wolf [****] in his oscillating 
flow study for tubes. 

With similarity to the pulsating flow, a lead phase angle for the 
pressure gradient over the velocity exists. Table 4.3.1 reveals the present 
lead phase angle as 68 degrees for the <o*=32.0 and Re ffiax =2000; this is 
compared to the 66 degrees found by Kurzweg [1985] at the same 
conditions. Table 4.3.1 also shows the lead phase angle results for the 
variable fluid property models. The graphical presentation of the pressure 
gradient variation for one oscillating cycle is illustrated in Figure 4.3.3, 
and sine wave pressure gradient curves were obtained. 
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Figure 4.3.1: 


Cross-sectional velocity profiles compared between present 
work with data from Kurzweg [1985] (fully-developed, 
incompressible oscillating flow, R 6^=2000, <■> =32.0) [OS-1]. 





Center-line pressure in the channel at different ot values 
(Incompressible oscillating flow, Re==2000, to =32.0) [OS-1]. 


Figure 4.3.2: 




Table 4.3.1: Comparison between present work and data from Kurzweg [1985] 

on the phase angle difference between cross-sectional average 
velocity and axial pressure gradient (fully-developed oscillating 
flow, R 6^=2000, 0=32.0) [OS-1 to OS-3]. 


«■ 

Phase Angle Difference (degrees) 

Kurzweg [1985] 

Present 

Incomp. 

Incomp. 

T.E. 

Comp. 

32.0 

66. 

68. 

65. 

66. 


Notes: 

Incomp. - Incompressible 
T.E. - Thermally Expandable 
Comp. - Compressible (M=.02) 








The examination of the mass flux indicated similar results as found 


in the pulsating flow analysis. For the incompressible model, the 
instantaneous mass flux is constant throughout the channel, (see Figure 
4.3.4a), while the compressible model indicated a delay of the input mass 
flux signal caused by the fluid density variation. Notice also that the 
existence of a mass wave front produced by the acceleration portions of 
the cycle oscillates through the channel (see Figure 4.3.4b). This wave 
front identifies the location of high density compression at any 
instantaneous <ot caused by the variable flow speed. 

The next parameter examined is the time-averaged apparent friction 
in the channel. Unlike the pulsating flow where the time-averaged value 
is the average value of a complete fluctuation cycle, the time-averaged 
value for oscillating flow is defined as the average value of the forward 
half of the cycle. As for the reverse half of the cycle, the time-averaged 
values are mirror images across the central channel of the forward half 
results. Figure 4.3.5 shows the time-averaged friction comparison between 
the different fluid property models plotted with the steady incompressible 
results as reference. Notice the close agreement between the steady 
incompressible results with the time-averaged oscillating friction 
predictions. This suggests that the oscillating flow is close to being 
quasi-steady. Table 4.3.2a lists the numerical friction values of Figure 
4.3.5, and substantial friction difference of about 50% is observed between 
the incompressible results and the thermally expandable results. The 
friction comparison between thermally expandable and compressible (M=.02) 
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Mass flux variation in the channel at <**=32.0 and Re lgz s2000: 

a) incompressible [OS-1]. 

b) compressible, M®.02 [OS-3]. 


Figure 4.3.4: 






models are close in value at the entrance, but diverges downstream at 15% 
difference. This frictional analysis expresses the importance to examine 
oscillating flow with the variable property models and reveals the limitation 
of the incompressible model. 

In Figure 4.3.6, the mean fluid temperature versus time is plotted for 
two axial locations: one at the central channel and one close to the 
entrance. The figure shows two oscillating temperature cycles for every 
oscillating velocity cycle which corresponds to other oscillating flow 
literature. The minimal difference in mean temperature is observed 
between different fluid property models. 

Figure 4.3.7 shows the time-averaged local heat flux in the channel 
predicted by the different fluid property models. With the wall/inlet 
temperature ratio set at 1.2, the heat flux predictions of the thermally 
expandable model are about 3% higher than the incompressible model. 
Identical time-averaged heat flux results were found between the thermally 
expandable model and the compressible model at M=.02 (see Table 4.3.2b). 

The same conclusion in the pulsating flow is observed for the 
oscillating flow. For variable property flows, the wall/inlet temperature 
ratio and the Mach number has a direct effect on the friction factor while 
a secondary effect on the heat transfer. 
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Figure 4.3.7: Comparison of the time-averaged local heat flux in the channel 

between different fluid property models (T^/t^I. 2, Re Ja =2000, 
o*=32.0) [0S-1, OS-2, OS-3]. 





Table 4.3.2: Comparison of the oscillating flow results between different fluid 

property models (Re^ZOOO, <i> =32.0, T^/r^l.2, Helium) [OS-1 to 
OS-3]: ^ 


a) time-averaged apparent friction. 


m 

Time-Averaged Local F* 

Diff. from T.E. (%) 

Incomp. 

mum 

Comp. 

Incomp. 

Comp. 

.01 

1.738 

3.283 

3.279 

-47.0 

-0.1 

.02 

1.409 

2.890 

2.860 

-51.2 

-1.0 

.04 

1.185 

2.523 

2.406 

-53.0 

-4.6 

.08 

1.037 

2.026 

1.725 

-48.8 

-14.9 


b) time-averaged heat flux. 


H 

Time-Averaged Local Q* 
(• 1000) 

Diff. from T.E. (%) 

Incomp. 


Comp. 

Incomp. 

Comp. 

.01 

.636 

.651 

.651 

-2.3 

0.0 

.02 

.487 

499 

.499 

-2.4 

0.0 

.04 

.370 

.381 

.381 

-2.9 

0.0 

.08 

.286 

.293 

.292 

-2.4 

-0.3 


Notes: 

Q* = -vc'oe/JY) ly_o 
F* - U * Re. / 24 
Incomp - Incompressible 
T.E. - Thermally Expandable 
Comp. - Compressible (M=.02) 
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Oscillating Flows fSPRE Boundary Conditions) : 

The two components under examination are the heater and the cooler 
of the NASA Space Power Research Engine (SPRE). In Appendix A-5, the 
run numbers OS-4 and OS-5 show the characteristic parameters of the 
cooler while the run numbers OS-6 and OS-7 represent the heater. Each 
component was solved under the incompressible fluid model and the 
compressible fluid model. 

Figure 4.3.8 shows the time-averaged local apparent friction of the 
heater and cooler compared with the results from the incompressible and 
compressible (M=.01) analysis. It can be seen that only a slight difference 
existed between the predictions of the two fluid property models. The 
same conclusion was found on the comparison of the time-averaged Nusselt 
number presented in Figure 4.3.9. 

The identical friction and heat transfer results between the 
incompressible and compressible models is caused by the characteristic 
parameters used. In the present analysis, the compressibility of the fluid 
is limited to the dependence of the temperature change and the flow speed. 
Since the wall/inlet temperature ratios are close to one and the Mach 
number is at a low .01, the conditions make these oscillating flows ideal 
incompressible property models. However, in the actual engine, the fluid 
compression and expansion caused by the out-of-phase movement between 
the displacer and the power piston will make an impact on the friction and 
the Nusselt number, which is not accounted for in this study. This also 
confirms the assumption Seume [1988] suggested that compressibility 
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initiated by the flow speed is insignificant compared to the volumetric 
change in the expansion and compression spaces. But to model the actual 
fluid compressibility caused by the volumetric change requires moving 
boundaries which is out of the current code's capacity range. 
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CHAPTER V 

CONCLUSIONS AND RECOMMENDATIONS 


The Partially-Parabolized Navier-Stokes equations have been used to 
analyze the compressible and unsteady, periodic internal flows. The 
computer code has been validated for steady flow conditions, and the 
present numerical predictions were in excellent agreement with data 
obtained from the Navier-Stokes equations and the boundary-layer 
equations reported in the literature. 

Upon examining the pulsating flow, the accuracy of the present 
numerical code was validated by the comparison with the closed form 
solutions from Siegel and Perlmutter [1962]. Excellent agreements were 
obtained for the incompressible fully-developed flow conditions and 
thermally developing with slug flow approximation. 

An interesting observation was made on the mass flux during the 
examination of the pulsating flow with variable fluid properties. The 
viscous damping of the fluid was diminishing the fluctuating component of 
the pulsating flow. Thus, if the channel is sufficiently long, the pulsating 
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flow will become steady-state downstream. This is an observation in which 


the incompressible model cannot obtain. 

Also found in the pulsating flow analysis is that the fluid 
compressibility has a stronger effect on the friction factor value than on 
the heat transfer coefficient (see Table 4.2.5 and 4.2.8). Furthermore, the 
difference between the steady state Nusselt number and the local time- 
averaged Nusselt number is less than 1.2% for the cases considered. This 
leads to the conclusion that pulsating flows contribute negligible heat 
transfer enhancements compared with the steady flows for the cases 
considered (see Appendix A-4). 

For the oscillating flow, the present fully-developed velocity profiles 
are compared with the analytical results from Kurzweg [1985] at <o*=32 and 
Re max =2000. Excellent agreement was obtained. From the variable property 
models where T v /t 0 =1.2, the substantial friction difference of about 50% was 
observed between the incompressible results and the thermally expandable 
results. This shows the importance to examine oscillating flow with the 
variable property models. 

For the SPRE conditions, the present code predicted similar friction 
and Nusselt number results from the incompressible and compressible 
models. However in the present analysis, the compressibility of the fluid 
is limited to the dependency of the temperature change and the flow speed, 
while neglecting the actual compression and expansion effects produced by 
the out-of-phase movement of the piston and the displacer. Further 
modification of the present code is needed to simulate the engine 
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conditions. Perhaps incorporating the moving inlet and exit boundaries 
and use of the full Navier-Stokes equations would predict results closer 
to the engine conditions. As an ultimate goal, it would be desirable to 
extend the present solution procedure to handle turbulent flows. 

Finally, two flow visualization videos were generated from the 
present solutions. The first video named "Pulsatile Flow Between Parallel 
Plates" animates the pulsating flows at three different Valensi numbers: .08, 
32.0, and 200.0. The fluid air travels at a mean Reynolds number of 2000 
while the amplitude of the velocity fluctuation is half of the mean velocity. 
In all cases, the fluid is heated by the constant wall temperature at the 
ratio of 1.2 to the inlet fluid temperature. 

The second video named "Oscillating Flow Between Parallel Plates" 
contains two parts of visualization animations. The first part shows three 
different Valensi numbers (.08, 32.0, and 100.0) with the operating maximum 
Reynolds number set at 2000. Helium is the selected fluid and is heated 
by the constant wall/inlet temperature ratio of 1.2. The second part 
contains the oscillating flow animations, at the Stirling engine cooler and 
heater conditions. The characteristic parameters are listed in Appendix A- 
5. 

Both pulsating and oscillating flows utilize a channel with a length 
of 70 times the height. Numerical calculations were done at every 5 
degrees, hence each cycle is animated by 72 frames of plots. The resulting 
parameters examined in both videos were the velocity profile, the 
temperature contour, the inlet velocity, the center-line pressure in the 


9 /jr 


channel, and the apparent friction and Nusselt number; each of these 
parameters is, animated with respect to time. 

The videos clearly reflect the quasi-steady behavior of the flow when 
the Valensi number is low (.08) but this behavior switches to non-quasi- 
steady as the Valensi number increases (32.0, 100.0 and 200.0). The videos 
also demonstrate the difference between incompressible, thermally 
expandable and compressible flows. 
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APPENDIX - A 


List of Flow Runs Examined 


Appendix A— 1 : Steady Flow Runs (Incompressible) 


Run 

Number 

Re. 

Fluid 
(Pr no.) 

Grid Size 
(x.y) 

Wall Temp, 9 

SI-1 

40 

Air (.72) 

50 X 91 

2.0 

SI-2 

100 

Air (.72) 

50 X 91 

2.0 

SI-3 

100 

Air (.72) 

50 X 91 

1.5 (upper) 
0.5 (lower) 

SI-4 

144 

Air (.72) 

60 X 40 

1.0 

SI-5 

1000 

Air (.72) 

60 X 40 

1.0 

SI-6 

2000 

Air (.72) 

107 X 29 

.67 

SI-7 

2000 

Air (.72) 

107 X 29 

1.0 

SI-8 

2000 

Air (.72) 

107 X 29 

1.5 


Appendix A-2: Steady Flow Runs (Thermally Expandable) 


Run 

Re. 

Fluid 

Grid Size 

Wall Temp, 6 

Number 

(Pr no.) 

(x,y) 

ST-i , 

144 

Air (.72) 

60 X 40 

0.5 

ST-2 

144 

Air (.72) 

60 X 40 

2.0 

ST-3 

144 

Air (.72) 

60 X 40 

5.0 

ST- 4 

2000 

Air (.72) 

107 X 29 

0.5 

ST- 5 

2000 

Air (.72) 

107 X 29 

0.67 

ST- 6 

2000 

Air (.72) 

107 X 29 

1.2 

ST- 7 

2000 

Air (.72) 

107 X 29 

1.4 

ST-8 

2000 

Air (.72) 

107 X 29 

1.5 

ST- 9 

2000 

Air (.72) 

107 X 29 

2.0 

ST- 10 

2000 

Air (.72) 

1 107 X 29 

5.0 


































































































Appendix A-3: Steady Flow Runs (Compressible) 


Run 

Number 

Re. 

Fluid 
(Pr no.) 

Grid Size 
(x,y) 

Wall Temp, 9 

Mach No. 

SC-1 

2000 

Air (1.0) 

107 X 29 

0.980 

0.25 

SC-2 

2000 

Air (1.0) 

107 X 29 

0.993 

0.25 

SC-3 

2000 

Air (1.0) 

107 X 29 

1.000 

0 25 

SC-4 

2000 

Air (1.0) 

107 X 29 

1.007 

0.25 

SC-5 

2000 

Air (.72) 

107 X 29 

0.67 

0.05 

SC-6 

2000 

Air (.72) 

107 X 29 

0.67 

0.10 

SC-? 

2000 

Air (.72) 

107 X 29 

0.67 

0.25 

SC-8 

2000 

Air (.72) 

107 X 29 

1.50 

0.05 

SC-9 

2000 

Air (.72) 

107 X 29 

1.50 

0.10 

SC-10 

2000 

Air (.72) 

107 X 29 

1.50 

0.25 



























































Appendix A-4: Pulsating Flow Runs 


Run 

Number 

cS 


T./T. 

Fluid 

Properties 

Mach no. 

PI-1 

0.08 

1.00 

2.0 

Incomp. 

— 

Pl-2 

8.00 

.778 

2.0 

Incomp. 

— 

Pl-3 

32.0 

.298 

2.0 

Incomp. 

— 

Pt-4 

200. 

.054 

2.0 

Incomp. 

— 

PI-5 

0.08 

1.00 

1.4 

Incomp. 

— 

pr-6 

32.0 

1.00 

1.4 

Incomp. 

— 

PI-7 

100. 

1.00 

1.4 

Incomp. 

— 

PT-1 

0.08 

1.00 

1.2 

T.E. 

— 

PT-2 

0.08 

1.00 

1.4 

T.E. 

— 

PT-3 

32.0 

1.00 

1.2 

T.E. 

— 

PT— 4 

32.0 

1.00 

1.4 

T.E. 

— 

PT-5 

1 100. 

1.00 

1.2. 

T.E. 

— 

PT-6 

100. 

1.00 

1.4 

T.E. 

— 

PC-1 

0.08 

1.00 

1.2 

Comp. 

0.05 

PC-2 

32.0 

1.00 

1.2 

Comp. 

0.05 

PC-3 

100. 

1.00 

1.2 

Comp. 

0.02 


Notes: All the above runs have the following parameters - 

a) Air as the fluid with Pr=72. 

b) x.y-grid is 107 X 29. respectively. 

c) Re_»2000. 

d) Time step ie 5 degrees. 
















































































Appendix A-5: Oscillating Flow Runs 


Run 

Number 

g/ 

Re*.. 

n 

T./t. 

Fluid 

Properties 

Mach no. 

OS-1 

32.0 

2.000 

.893 

1.2 

Incompressible 

-- 

OS-2 

32.0 

2,000 

.893 

1.2 

Thermally Expandable 

-- 

OS-3 

32.0 

2,000 

.893 

1.2 

Compressible 

.02 

OS-4 

350.0 

30.000 

1.22 

0.956 

Incompressible 

-- * 

OS-5 

350.0 

30,000 

1.22 

0.956 | 

Compressible 

.01 

OS-6 

88.0 

16.500 

2.68 

1.048 

Incompressible 

-- 

OS-7 

88.0 

16.500 

2.68 

1.048 

Compressible 

' .01 


Notes: All the above runs have the following parameters - 

a) Helium as the fluid with Pr*.72. 

b) x,y-grid is 107 X 29, respectively. 

c) Channel length is L/H=70.0. 

d) Time step is 5 degrees. 
















































APPENDIX - B 

Explanation of Input Parameters for the Computer Code 
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This variable list contains the input parameters used by the present 
computer code and their explanations. The sequence of these variables 
are in the same order as they appear in the READ statement. 


HIGHT 

THIGHT 

HSTEP 

XSTEP 

FLARE 

DMLFCT 

LPNS 

KPNS 

JPNS 

NLMT 

LPOP 

NJ 

INV 


Channel inlet height. 

Total channel height including the wall's thickness. If 
the wall's thickness is not considered in the analysis, 
set equal to HIGHT. 

Step height of the sudden expansion. For Straight 
Parallel Plates Channel(SPPC), set equal to 0.0 . 

Streamwise location of the step for SPPC, set equal to 

0.0 . 


Constant for FLARE approximation. Set equal to 0.0 . 

Constant used for the set up of the x-grid. Normally, 
set equal to 1.01 - 1.20 . 

Designates the axial station index (MCOUNT) after 
which the governing equations (PPNS or NS) are to be 
solved. For most cases, set equal to 1 . 

Designates MCOUNT after which the governing 
equations need not be solved. 

Designates MCOUNT after which Global Iteration in 
Space(GIS) is desired. In general, set to 2 . 

Safety parameter designating the maximum allowable 
MCOUNT, set equal to KPNS. 

For external pressure gradient flows, this represents 
the number of freestream u-velocity input. For 
internal fluid flows, set to 0 . 

Total number of grid points in the y-direction inside 
the channel. 

For external pressure gradient flows, this represents 
the number of freestream v-velocity inputs. Set to 0 
for internal flows or when no v-velocities are 
prescribed. 
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INT This represents the number of wall temperature inputs 

on the upper or lower boundaries. Set to 0 for 
isothermal or constant heat flux boundary conditions. 

ZAP1 Output selections: For abbreviated output, set equal to 

1 ; for detailed output at each streamwise station, set 
to -1 . 

GLOBE Beginning cycle number for this run. Set equal to 

min. 

TOLERC No longer used. 

Tolerance on the streamwise pressure gradient to be 
used in secant procedure when estimating the initial 
pressure field for the first cycle. 

XE Safety parameter designating the maximum axial 

distance in length beyond which no calculation will be 
done. 

PR Prandtl number. 

CPS Specific heat. 

TEST Value used to check for edge of B.L., typically 0.9995 

for fully-developed flow and 0.995 for developing flow. 

US The u-velocity used for nondimensionalization. 

XMUS The absolute viscosity at the channel's inlet used for 

non dimensionalization. 

RHOS The fluid density at the channel's inlet used for 

nondimensionalization. 

VW Normal velocity at the wall. Set equal to 0.0 for no- 

slip conditions at the wall. 

UREF Freestream u-velocity that will vary with axial 

distance, however, US remains fixed. Set equal to the 

freestream value at first axial station. 

TWL Wall temperature for the lower isothermal boundary. 

TWU Wall temperature for the upper isothermal boundary. 

For the symmetrical analysis, set equal to any value. 
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MIN 

MAX 

IUPDAT 


NPOINT 


LTEMP 


MP1 

MPC(J) 

NP1 

NWALL 


JSTEP 


MCSTEP 


Beginning GIS. In general, set to 1 . 

Final GIS. 

Set to 1 if the best known values are used in 
evaluating convective coefficients for the pressure 
gradients, FI and F2. Otherwise set equal to 0 . In 
general, set equal to 0 . 

The number of y-grid points used inside of each 
channel wall. If the wall thickness is not considered 
in the analysis, set equal to 0 . 

Number of x-stations from the inlet for the wall 
temperature to reach the desired constant wall 
temperature. For constant wall temperature starting 
at the inlet, set equal to 0 . 

The number of MCOUNTs to be read in to determine 
detailed printout stations. 

J=1,MP1 . The MCOUNTs of streamwise stations where 
detailed printout is desired. 

Number of x distance from the inlet to be read in to 
determine detailed printout locations. 

For asymmetric test cases: 

Set to 2 if one-sided differencing is to be used to 
evaluate FI and F2 near wall. 

Set to -2 if regular differencing is to be used. 

For other geometries (including SPPC): 

Set to 1 for one-sided differencing. 

Set to -1 for regular differencing. 

Grid-point index in y-direction for the first point 
below the sudden expansion step. For the SPPC 
analysis, set equal to 1 . 

Station index (MCOUNT) for the first station 
downstream of the sudden expansion step. For the 
SPPC analysis, set to 1 . 

Set to 1, if the energy equation is to be solved, 
otherwise set to 0 . For variable fluid properties, set 
NTEMP equal to 1 . 
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NTEMP 


NWALLT 


NEXTRP 

XP3(J) 

NPRINT 


LS0R1 


LSOR2 


LEAD 


FAC 


Set to 0, if both the upper and lower boundaries are 
isothermal. 

Set to 1, if the temperatures on both boundaries are 
not constant. 

Set to >1, if only the lower boundary is NOT 
isothermal. 

Set to <0, if only the upper boundary is NOT 
isothermal. 

Set to 1, if extrapolation is to be used to evaluate FI 
and F2 near wall, otherwise set to 0 . 

J=1,NP1 . The x-distances from the inlet of the 
channel where detailed printout is desired. 

Set to 0 if abbreviated output is required when 
solving the Poisson equation for pressure. Set to 1 
for detailed printout. ZAP, if set equal to 1, 
overrides NPRINT and no details are printed. 

For the steady flow analysis: 

If set to 1, calculation of the beginning GIS (MIN) 
starts with solving the Poisson equation. In this 
case, the marching sweep for this GIS must have been 
already calculated and the results stored on unit 9. 
Choice of LSOR1 is related to value of LSOR. Set to 0 
for first GIS. 

For the unsteady flow analysis: 

Set to 0 for the first time cycle. 

If set equal to 1, the final GIS (MAX) completes with 
the calculation of the Poisson equation for pressure. 

If set equal to 0, the final GIS finishes with the 
marching integration sweep calculation but the Poisson 
equation is not solved. 

This variable specifies the singularity condition of the 
u-velocity at the point just ahead of the leading edge 
of the plate. When set to 0, this u-velocity is zero. 
When set to 1, this u-velocity is set to UREF. In the 
present study (SPPC), LEAD is set to 0. 

The over-relaxation factor used on the pressure 
update after solving the Poisson equation by the 
method of SOR by points. Values used range from 1.5 
to 1.97 . 
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TOL 


XKWL 

FROMINLT 


RC 

ERRMAX 

NIRROT 

IRROT 


NDOWN 


NOBLK 


NBLK 


No longer used. 

When positive, this represents tolerance on the total 
mass flow rate when making block adjustments on 
pressure. This option used only after convergence 
has become montonic. If block adjustments are not 
required, set to any negative value. Make sure TOL 
is small when positive to avoid imposing oscillations. 

The thermo-conductivity of the walls If the 

analysis does not include the walls, set equal to PR. 

This variable is for print-out purpose only. Set to 
the x-value where the current results are started 
from. In general set this equal to 0.0 . 

Critical mesh Reynolds number. Set to 1.9 . 

Mass convergence criteria. Set to around .10 . 

Set to 2 . 

Set to 0 if the freestream u, v-velocities are specified. 
Set to 1 if only freestream u values are specified and 
an irrotational outer-edge boundary condition is to be 
used on the v-velocities. The second option does not 
work well. 

Specifies the downstream pressure boundary condition. 
If set to: 

0 - no downstream boundary condition is imposed on 

the pressure. 

1 - constant axial pressure gradient condition is 

imposed. 

-1 - constant normal pressure gradient condition is 

imposed. 

When set to 0, the pressure block adjustment is used 
in all the GIS. If set to 1, the pressure block 
adjustment is used in the GIS up to the iteration 
specified in NBLK, after that no pressure block 
adjustment is used. For the variable properties 
analysis, set NOBLK to 1 . 

The number of GIS that the pressure block adjustment 
is used. NBLK is ignored when NOBLK is set to 0 . 
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NSUTH 


This variable specifies the method used to determine 
the viscosity of the fluid. If set to 1, the 
Sutherland's Law is used. If set to 0, the Power Law 
is used. For the constant properties case, set NSUTH 
to 0 . 

NCOMP For compressible flow, set equal to 1 . 

For thermally expandable or incompressible flow, set 
equal to 0 . 

XBEGIN This variable represents the x coordinate value at 

which the calculation begins. 

SUTH The Sutherland's constant in degree Kelvin used to 

determine the value of the fluid viscosity. (Example: 
for air SUTH=110.4 K) SUTH is ignored when NSUTH 
is set to 0 . 

ALPHA The exponent a in the Power Law is used to determine 

the fluid viscosity. Set ALPHA to 0.0 for 
incompressible analysis. ALPHA is ignored when 
NSUTH is set to 1 . 

BETA The exponent p in the Power Law is used to determine 

the fluid conductivity. Set BETA equal to 0.0 for 
incompressible analysis. BETA is ignored when NSUTH 
is set equal to 1 . 

GAMMA The exponent y in the perfect gas relationship is used 

to determine the fluid density. Set GAMMA equal to 
1.0 for variable density (thermally expandable or fully 
compressible). Set GAMMA equal to 0.0 for 
incompressible analysis. 

XMACH The reference mach number of the flow. This variable 

is ignored when NCOMP is equal to 0 . 

CPRATIO The ratio of specific heats of the fluid. This variable 

is ignored when NCOMP is equal to 0 . 

NSTEAD This variable specifies the type of flow to be 

considered. Set to 1 for steady flow. Set to 0 for 
unsteady flow (pulsatile or oscillating). 

N2DPRT For print-out of 2-D velocity and fluid properties, set 

to 1.0 . 




UB( J) 
VB(J) 

NCYCLE 

NDEGRE 

NSLUG 

NOSCIL 

NSKIP 

VALENSI 

UMEAN 


Read only when the flow is steady (NSTEAD=1). 
J=1,NJ+1. The inlet u-velocity profile at the upstream 
boundary. 

Read only when the flow is steady (NSTEAD=1). 
J=1,NJ+1. The inlet v-velocity profile at the upstream 
boundary. Note that the axial location of this profile 
is slightly different than that for the u-velocity 
profile due to the staggered grid format. 

Read only when the flow is unsteady (NSTEAD=0). 

The number of pulsating (or oscillating) cycles to be 
made in this run. 

Read only when the flow is unsteady (NSTEAD=0). 

The degree increment used in the calculation (One 
cycle has 360 degrees). Set NDEGRE greater than or 
equal to 5 . 

Read only when the flow is unsteady (NSTEAD=0). 

Set to 1 for the slug flow approximation when 
determining the temperature field. Set to 0 for 
simultaneously developing flow. The slug flow 
approximation works only for the incompressible flow. 
Set to 0 when making a variable properties analysis. 

Set to 1 for oscillating flow. Set to 0 for pulsating 
flow. 

Number of time step to be skipped around the zero 
velocity for the oscillating flow. Input even number 
only. 

Read only when the flow is unsteady (NSTEAD=0). 
Valensi number. 

Read only when the flow is unsteady (NSTEAD=0). The 
time average value of the inlet u-velocity. For 
Oscillating flow, set UMEAN equal to 0.0 . Note that for 
unsteady flow the inlet velocities are assumed uniform 
in profile. 
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DELTU 


MSEP1 

MSEP2 


MCDOWN 
NS FULL 

NHYBRD 


Read only when the flow is unsteady (NSTEAD=0). The 
amplitude of the fluctuating u-velocity at the inlet of 
the channel. (Example: If UMEAN & DELTU are set to 

1.0 & .5, respectively, the inlet u-velocity will vary 
from .5 to 1.5 sinusoidally through the cycle). 

Set to the MCOUNT value beyond which the 2-D 
u,v-velocity arrays are stored. For separated, 
incompressible, and steady flow, set equal to the 
MCOUNT value slightly ahead of the separation point. 
For steady incompressible flow in SPPC, set equal to 
KPNS. For other flows (variable properties, pulsating, 
or oscillating) set equal to 1 . 

Set to the MCOUNT value beyond which the 2-D 
u,v-velocity arrays are not stored. For separated, 
incompressible, and steady flow, set equal to the 
MCOUNT value slightly beyond the reattachment. For 
other flows (steady incompressible flow in SPPC, 
variable properties, pulsating, or oscillating) set equal 
to KPNS. For full Navier-Stokes calculation, the 
specified values of MSEP1 and MSEP2 are overridden 
within the program (no changes were made in the full 
N.S. equations to run variable properties or unsteady 
flows). 

The MCOUNT value at the downstream. Set equal to 
KPNS. 

Set equal to 0 for solving the Partially Parabolized 
Navier-Stokes (PPNS) equations. Set to 1 if the full 
Navier-Stokes (NS) equations are to be solved. Note 
that even when full NS equations are to be solved, 
NSFULL must be set to 0 for the first GIS. No 
modifications were made to update the full NS 
equations for unsteady, thermally expandable, or 
compressible flow analysis. For the present study, 
set to 0 . 

Set equal to 1 to use the hybrid differencing scheme. 

If set to 0, pure upwind or central differencing is 
used depending on the value of the mesh Reynolds 
number. In general, set to 1 . 


& /St 


NSTEP 

MSEP11 

MSEP22 

XU(J) 

YU(J) 

XV(J) 

YV(J) 

MINT 

MAXT 

NENGY 

IBL 


Set equal to 0, when analyzing the external or 
internal asymmetric channel flow. Set to 1 for 
internal symmetric channel flow. 

Set to the MCOUNT value where beyond which the 2D 
u,v-velocities were stored from the previous GIS. In 
general, set equal to MSEP1. 

Set to the MCOUNT value where beyond which the 2D 
u,v-velocities were Not stored from the previous GIS. 
In general, set equal to MSEP2. 

Not required for internal flows. 

J=1,LP0P. The x distances for freestream u-velocity 
input. (Note: No attempts were made to run an 

external flow analysis under the unsteady, thermally 
expandable, or compressible conditions). 

Not required for internal flows. 

J=1,LP0P. The freestream u-velocity values 
corresponding to XU(J) values. 

Not required for internal flows or if INV is set to 0 . 
J=1,INV. Represents the x distance for freestream 
v-velocity input. 

Not required for internal flows or if INV is set to 0 . 
J=1,INV. Represents the freestream v-velocity values 
corresponding to XV(J) values. 

Beginning GIS number for solving the energy 
equation. 

Final GIS number for solving the energy equation. 

Set to 1 if only the energy equation is to be solved 
(this option requires that the u,v-velocities have been 
solved in previous run); otherwise set to 0 . 

Set equal to 1, if the Boundary Layer (B.L.) or the 
compressible energy equation is used, otherwise set 
to 0 . For Thermally expandable or Compressible flow 
analysis, set IBL equal to 1 . 




NQWL 

NQWH 

QWL 

QWH 

TREF 

TT(1,J) 

XTWL(J) 

YTWL(J) 

XTWU(J) 

YTWU(J) 

NG3(I) 


Set to 1 for constant heat flux at the lower wall, 
(currently, this works only for incompressible flow), 
otherwise set equal to 0 . 

Set to 1 for constant heat flux at the upper wall, 
(currently, this works only for incompressible flow), 
otherwise set equal to 0 . 

This is the constant heat flux value at the lower wall, 
(currently, this works only for incompressible flow), 
if NQWL is set to 0, then QWL is ignored. 

This is the constant heat flux value at the upper wall, 
(currently, this works only for incompressible flow), 
if NQWH is set to 0, then QWH is ignored. 

The reference stagnation inlet temperature (in Kelvin). 

J=1,NJ. The inlet temperature (in Kelvin) profile at 
the upstream boundary. 

Not required for Isothermal or constant heat flux 
boundary. J=1,INT. The x distances for lower 
boundary temperature inputs. 

Not required for Isothermal or constant heat flux 
boundary. J=1,INT. Lower boundary temperature 
values corresponding to the x distances, XTWL(J), 
values. 

Not required for Isothermal or constant heat flux 
boundary. J=1,INT. The x distances for upper 
boundary temperature inputs. 

Not required for Isothermal or constant heat flux 
boundary. J=1,INT. Upper boundary temperature 
values corresponding to the x distances, XTWU(J), 
values. 

I=MIN,MAX. The number of Gauss-Siedel sweeps to be 
carried out when solving the Poisson equation for 
pressure. For the first GIS, set to a small number 
(around 5). For later cycles, increase this number 
gradually. 




NT3(I) 


No longer used. 

I=MIN,MAX. The number of times the pressure is 
revised at each axial station during each conventional 
iteration of the SOR method when solving the Poisson 
equation for pressure. 

FAC13(I) I=MIN,MAX. The under-relaxation factor to be used on 

the pressure gradients after each marching 
integration sweep. Start with a small fraction and 
increase gradually but not to excess one. 

FAC23(I) I=MIN,MAX. The relaxation parameter for velocity 

corrections; 1.0 was used in this study. 

FAC33(I) I=MIN,MAX. The relaxation parameter for fluid 

properties corrections; 1.0 was used in this study. 

*** Note: One data set of NG3, NT3, FAC13, FAC23, & FAC33 values 

is required for each GIS. 
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